2
$\begingroup$

I would like to solve the following equations simultaneously and numerically for all $X, Y, Z, W$ where i = 1:Nw, j = 1:Nl, k = 1:K.

$W_\text{net1}$, $W_\text{net2}$, $m_\text{net1}$, and $m_\text{net2}$ are given.

How can I apply solve() or any other function to numerically find the values of $X, Y, Z, W$?

enter image description here

Here is the code for implementing the above equation in matlab.

Nw = 2; Nl=3; K=2;
ivec = 1 : Nw;
jvec = 1 : Nl;
W_net1 = randi([-123, 123], Nw, K);
W_net2 = randi([-123, 123], Nl, K); 
m_net1 = randi([-123, 123], Nw, K); 
m_net2 = randi([-123, 123], Nl, K); 

for k = 1 : K
    for i = ivec
        X(i, k) = 2*(1 - 2*W(i,k))/((1 - 2*W(i,k))*(1 + W_net1(i,k)) + W(i,k)*W_net1(i,k)*(1 - (2*W(i,k))^m_net1(i,k)));
        ii = setdiff(ivec, i);
        tW1 = prod( 1 - X(ii, k) );
        tW2 = prod( 1 - Y(jvec, k) );
        W(i,k) = 1 - tW1 * tW2;
    end
    for j = jvec
        Y(j, k) = 2*(1 - 2*Z(j,k))/((1 - 2*Z(j,k))*(1 + W_net2(j,k)) + Z(j,k)*W_net2(j,k)*(1 - (2*Z(j,k))^m_net2(j,k)));
        i = ivec;
        tZ1 = prod(1 - X(i, k));
        jj = setdiff(jvec, j);
        tZ2 = prod(1 - Y(jj, k));
        tZ3 = tZ2 * tZ1;
        Z(j,k) = 1 - tZ3;
        S(j,k) = Y(j, k) * tZ3;
    end
end
Improve this question
$\endgroup$

1 Answer 1

Reset to default
3
$\begingroup$

One way of approaching the problem you posted is to use fsolve.

You can start by creating a file with a function that returns the difference $dX$, $dY$, $dW$ and $dZ$ between an estimate of $X$, $Y$, $W$ and $Z$ and the respective expressions for these variables, $f_X$, $f_Y$, $f_W$ and $f_Z$ that you presented in your post: $$\begin{align} dX(i,k) &= f_X(X,Y,W,Z) - X(i,k) \\ dY(j,k) &= f_Y(X,Y,W,Z) - Y(j,k) \\ dW(i,k) &= f_W(X,Y,W,Z) - W(i,k) \\ dZ(j,k) &= f_Z(X,Y,W,Z) - Z(j,k) \end{align}$$

However, in order to call this function it's easier if you map the variables into a matrix $U$ of size $Nw + Nl + Nw + Nl$ by $K$. This way, you can use the matrix $U$ as the function argument. The parameters can be passed to the function by using function handles.

After creating the file with the function, your main program should define an initial guess (can it be all zeros, for instance?) and call fsolve.

Edit:

In order to define the function, you'll need to create a file called fun_dif.m with the following content:

function dU = fun_dif(U, Nw, Nl, K, W_net1, W_net2, m_net1, m_net2)

ivec = 1 : Nw;
jvec = 1 : Nl;

X(ivec, 1:K) = U(ivec, 1:K);
Y(jvec, 1:K) = U(Nw + jvec, 1:K);
W(ivec, 1:K) = U(Nl + Nw + ivec, 1:K);
Z(jvec, 1:K) = U(Nw + Nl + Nw + jvec, 1:K);
dU = zeros(Nw + Nl + Nw + Nl, K); % initialize the difference vector

for k = 1 : K
    for i = ivec
        dU(i, k) = 2*(1 - 2*W(i,k))/((1 - 2*W(i,k))*(1 + W_net1(i,k)) + W(i,k)*W_net1(i,k)*(1 - (2*W(i,k))^m_net1(i,k))) - X(i, k);
        ii = setdiff(ivec, i);
        tW1 = prod( 1 - X(ii, k) );
        tW2 = prod( 1 - Y(jvec, k) );
        dU(Nl + Nw + i, k) = 1 - tW1 * tW2 - W(i,k);
    end
    for j = jvec
        dU(Nw + j, k) = 2*(1 - 2*Z(j,k))/((1 - 2*Z(j,k))*(1 + W_net2(j,k)) + Z(j,k)*W_net2(j,k)*(1 - (2*Z(j,k))^m_net2(j,k))) -  Y(j, k);
        i = ivec;
        tZ1 = prod(1 - X(i, k));
        jj = setdiff(jvec, j);
        tZ2 = prod(1 - Y(jj, k));
        tZ3 = tZ2 * tZ1;
        dU(Nw + Nl + Nw + j, k) = 1 - tZ3 -  Z(j,k);
        S(j,k) = Y(j, k) * tZ3;
    end
end

Note that it's mostly your code with the computation of the difference $dU$ between the estimate from the formulas and the current value $U$.

Then, in your main program you can have

Nw = 2;
Nl=3;
K=2;
W_net1 = randi([-123, 123], Nw, K);
W_net2 = randi([-123, 123], Nl, K); 
m_net1 = randi([-123, 123], Nw, K); 
m_net2 = randi([-123, 123], Nl, K); 

f = @(U)fun_dif(U, Nw, Nl, K, W_net1, W_net2, m_net1, m_net2); % create a function handle
U0 = 0.5 + zeros(Nw + Nl + Nw + Nl, K); % initial estimate is 0.5 for all matrices
UU = fsolve(f,U0); % UU will have the final solution
X(ivec, 1:K) = UU(ivec, 1:K);
Y(jvec, 1:K) = UU(Nw + jvec, 1:K);
W(ivec, 1:K) = UU(Nl + Nw + ivec, 1:K);
Z(jvec, 1:K) = UU(Nw + Nl + Nw + jvec, 1:K);

Edit:

It might be worthy to note that $$ Z(j,k) = 1 - \frac{S(j,k)}{Y(j,k)}$$ and $$ W(i,k) = 1 -\frac{S(j,k) (1 - Y(j,k))}{Y(j,k) (1 - X(i,k))} = 1 - \bigg(1 - Z(j,k) \bigg) \frac{1 - Y(j,k)}{1 - X(i,k)}$$

Improve this answer
$\endgroup$
9
  • $\begingroup$ Thank you so much for your reply. However, I am not sure if I got you quite well. Could you please tell me how I can define f here? Thank you very much in advance! $\endgroup$
    – Susan
    May 13, 2019 at 14:57
  • $\begingroup$ @Susan I haven't tested the code, so it might have bugs. Nevertheless, I tried to present the idea of the program. Please let me know if you have any questions. $\endgroup$
    – Ertxiem - reinstate Monica
    May 14, 2019 at 1:39
  • $\begingroup$ Thank you very much!!! The code works well and now I understand the approach you took to solve the problem. I DO appreciate your help . $\endgroup$
    – Susan
    May 14, 2019 at 15:22
  • $\begingroup$ You're welcome. I'm glad I could help. $\endgroup$
    – Ertxiem - reinstate Monica
    May 14, 2019 at 20:30
  • $\begingroup$ I would do a full numerical approach like the one I did above. Following your comment above, the 'W_net2`are the "variables", together with $X$, $W$, $Y$ and $Z$; and $S$ are treated as fixed parameters. I'm sure you can adapt my code to work like this. Perhaps there is a simpler way to solve this if (some of) the equations above have analytical solutions, but I'm unable to see how to get them. $\endgroup$
    – Ertxiem - reinstate Monica
    May 16, 2019 at 0:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.