# Solving $n$ coupled equations numerically in Matlab

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$$?

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


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)}$$

• 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! – Susan May 13 '19 at 14:57
• @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. – Ertxiem - reinstate Monica May 14 '19 at 1:39
• 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 . – Susan May 14 '19 at 15:22
• You're welcome. I'm glad I could help. – Ertxiem - reinstate Monica May 14 '19 at 20:30
• 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. – Ertxiem - reinstate Monica May 16 '19 at 0:55