Solving system of n equations

Hi MATLAb guys,

I am stucking at some point and need your help.

Do you know how it is possible to find a relationship between two variables in one equations? for example y = ax+b, (a,b) are given, (x,y) are variables. I would like matlab give me x = (y-b)/a.

I have wrote the following code, but it doesn't work out well. Any idea how to fix it? Thanks in advance

Input: n_W, n_L, W_net1, m_net1, m_net2, beta

Output: W_net2

funX=@(W,   W_net1,   m_net1) (2*(1 - 2*W))/((1 - 2*W)*(1 + W_net1) + W*W_net1*(1 - (2*W)^m_net1)); 
funY=@(Z,   W_net2,   m_net2) (2*(1 - 2*Z))/((1 - 2*Z)*(1 + W_net2) + Z*W_net2*(1 - (2*Z)^m_net2));     
funW=@(X,Y,n_W,n_L) (1 - ((1 - X)^(n_W - 1))*(1 - Y)^n_L);                                           
funZ=@(X,Y,n_W,n_L)(1 - ((1 - X)^n_W)*(1 - Y)^(n_L - 1));                                             
funS=@(X,Y,n_W,n_L,beta)(Y*((1 - Y)^(n_L - 1))*(1 - X)^(n_W))-beta;
fun1=@(X,Y,W,Z,n_W,n_L,W_net2, m_net2, W_net1, m_net1, beta) ([funX(W,W_net1,m_net1); funY(Z,W_net2,m_net2) ; funW(X,Y,n_W,n_L); funZ(X,Y,n_W,n_L) ; funS(X,Y,n_W,n_L,beta)]); 

n_W = 3;
n_L = 4;
beta = 0.7;
m_net2 = 4;
W_net1 = 8;
m_net1= 4;

fun2=@(P) (P-fun1(P(1),P(2),P(3),P(4),n_W,n_L,beta,m_net1,W_net1,m_net2));
InitialGuess=[0;0;0;0];
fsolve(fun2,InitialGuess)

P.S. The approach is

1) knowing beta, from S ===> relationship between X and Y can be found.

2) from X and W ===> X and Y can be found numerically

3) based on Y, W_net2 can be derived

[Walter Roberson]

In fun2, you call fun1() with fewer parameters (10) than fun1 needs (11), and the parameters you provide are in a different order than fun1 uses. The missing parameter is W_net2 .

If you provide a W_net2 and you re-order the parameters in the call to match what are expected, then fun1 produces a result which is 5 x 1. In fun2 you subtract that from P, but P is only 4 x 1.

[Susan]

Thank you so much for your reply.

W_net2 is an output that I am looking for. Now, with the modification that you mentioned, I have got the following error:

Undefined function or variable 'W_net2'.

Any idea?

Thanks again

[Walter Roberson]

My investigations seem to be indicating that when W_net2 is

-1/1048576/((-(p3-1)^3*(64*p3^4+32*p3^3+16*p3^2+8*p3+7)^6*(64*p3^4+32*p3^3+16*p3^2+8*p3+9)^6)^(1/2)*(64*p3^4+32*p3^3+16*p3^2+8*p3+7)^3)^(1/3)*(p3^4+1/2*p3^3+1/4*p3^2+1/8*p3+9/64)^9*((343/2+1/2*(-(p3-1)^3*(64*p3^4+32*p3^3+16*p3^2+8*p3+7)^6*(64*p3^4+32*p3^3+16*p3^2+8*p3+9)^6)^(1/4)+131072*p3^12+196608*p3^11+196608*p3^10+163840*p3^9+141312*p3^8+92160*p3^7+52736*p3^6+27648*p3^5+14304*p3^4+5296*p3^3+1848*p3^2+588*p3)*((-(p3-1)^3*(64*p3^4+32*p3^3+16*p3^2+8*p3+7)^6*(64*p3^4+32*p3^3+16*p3^2+8*p3+9)^6)^(1/2)*(64*p3^4+32*p3^3+16*p3^2+8*p3+7)^3)^(1/3)+262144*(p3^4+1/2*p3^3+1/4*p3^2+1/8*p3+7/64)^3*(-(p3-1)^3*(64*p3^4+32*p3^3+16*p3^2+8*p3+7)^6*(64*p3^4+32*p3^3+16*p3^2+8*p3+9)^6)^(1/4)+((-(p3-1)^3*(64*p3^4+32*p3^3+16*p3^2+8*p3+7)^6*(64*p3^4+32*p3^3+16*p3^2+8*p3+9)^6)^(1/2)*(64*p3^4+32*p3^3+16*p3^2+8*p3+7)^3)^(2/3))/(9/36028797018963968*(-(p3-1)^3*(64*p3^4+32*p3^3+16*p3^2+8*p3+7)^6*(64*p3^4+32*p3^3+16*p3^2+8*p3+9)^6)^(3/4)+(p3^4+1/2*p3^3+1/4*p3^2+1/8*p3+9/64)^3*(-6122493/68719476736*p3+p3^27-3/4*p3^25-26671987/67108864*p3^9-12826905/67108864*p3^8-45213609/536870912*p3^7-36988441/1073741824*p3^6-25785699/68719476736*p3^2-94499583/68719476736*p3^3-28045983/2147483648*p3^5-11/4*p3^24-195/32*p3^23-693/64*p3^22-1791/128*p3^21-4191/256*p3^20-70857/4096*p3^19-68035/4096*p3^18-117387/8192*p3^17-47661/4096*p3^16-143977/16384*p3^15-809247/131072*p3^14-132135/32768*p3^13-656277/262144*p3^12-24320193/16777216*p3^11-6581607/8388608*p3^10-19832313/4294967296*p3^4-31/274877906944*(-(p3-1)^3*(64*p3^4+32*p3^3+16*p3^2+8*p3+7)^6*(64*p3^4+32*p3^3+16*p3^2+8*p3+9)^6)^(1/2)+49/2097152*(p3^4+1/2*p3^3+1/4*p3^2+1/8*p3+9/64)^3*(-(p3-1)^3*(64*p3^4+32*p3^3+16*p3^2+8*p3+7)^6*(64*p3^4+32*p3^3+16*p3^2+8*p3+9)^6)^(1/4)-1180531/68719476736))/(p3^12+3/2*p3^11+3/2*p3^10+5/4*p3^9+69/64*p3^8+45/64*p3^7+103/256*p3^6+27/128*p3^5+447/4096*p3^4+331/8192*p3^3+231/16384*p3^2+147/32768*p3-1/262144*(-(p3-1)^3*(64*p3^4+32*p3^3+16*p3^2+8*p3+7)^6*(64*p3^4+32*p3^3+16*p3^2+8*p3+9)^6)^(1/4)+343/262144)

then the equations can hold, provided that p3 <= 1 . That is, in at least one combination of choices of solutions, the equations seem to have a redundancy . However, when I put it through automatic solvers, the solvers tell me that there is no solution. That potentially hints that if enough back-substitution was done, that there would be a division by 0 on all of the paths.

[Susan]

Thanks again for your reply. I appreciate your time and help.

I see.

Do you think it is possible to follow this approach and get any result for Y or W_net2?

1) knowing beta, from S ===> relationship between X and Y can be found.

2) from X and W ===> X and Y can be found numerically

3) based on Y, W_net2 can be derived

If so, how can I implement this approach in MATLAB? Thanks.

[Susan]

Hi Walter,

Just a quik question for you. p3 = P(3)? Actually all of P(1)-P(4) should be <=1 and >=0. However when I solve the equations, P(1) is negative and P(2)-P(4) >=1.

Any idea?

[Walter Roberson]

After the adjustment of your code in fun2 to add a W_net2 parameter to the call to fun1, and re-order the parameters in the call to match what fun1 expects, then when you call fun2 you effectively have five unknowns, p(1), p(2), p(3), p(4), and W_net2, and calling fun2 returns five expressions. Therefore you can treat fun2 as generating five equations in five unknowns, and try to solve (closed form) or vpasolve() (symbolic numbers) or fsolve() (numeric)

When I ask to solve() or vpasolve(), MATLAB tells me that there is no solution.

When I move the problem over to Maple and proceed with step by step elimination of variables, the expressions can get rather complicated, and multiple choices can arise. If I isolate p1, p2, p4, W_net2 in that order, then I am left with an expression in p3 that needs to be solved. The expression involves the roots of negative numbers when p3 > 1. When I plot below that, the graph seems to show thousands of roots of low magnitude, down at the numeric noise level. When I examine any particular p3 value more carefully, I am told that the exact value of the expression is 0. This all suggests that at least under some choices of p1, p2, p4, W_net2, that one of the equations is redundant. One of the solutions might be near p1 = 0.00620298344073486, p2 = -0.273330608872044, p3 = -1.59633120075406, p4 = -1.02635998232398, W_net2 = -1.25386620069808

[Walter Roberson]

Yes, p3 = P(3) . It is easier to work symbolically with named variables rather than indexed variables.

[Walter Roberson]

My tests are pretty sure that there is no solution when p1 to p4 are constrained to 0 to 1. However, in those circumstances I only tested with W_net2 in the range -1000 to +1000; the tests suggest that the best W_net2 in that range is about -0.787 . There is a fairly noticeable sum-of-squares residue of about 0.17 for all p1 to p4 constrained to 0 to 1. By way of contrast, near the positions I indicated above, the sum-of-squares residue is less than 2E-14

[Susan]

@Walter: Thank you VERY much for your thorough reply. It is awesome to see how volunteers like you spend their time to help others. I DO appreciate your time.

I have got another question. It would be greatly appreciated if you could kindly answer that. I am pretty new to this area and sorry in advance if my questions are silly.

1) As you mentioned, when I call fun2 I effectively have five unknowns, p(1), p(2), p(3), p(4), and W_net2, and calling fun2 returns five expressions. Is it anyway to access these five expression? when I try, I get

fun2 =

function_handle with value:

@(P)(P-fun1(P(1),P(2),P(3),P(4),P(5),m_LAA,n_L,Wmin_WiFi,m_WiFi,n_W,betaLAA))

2) when I modified the functions as follows

funX=@(W,W_net1,m_net1) (2*(1 - 2*W))/((1 - 2*W)*(1 + W_net1) + W*W_net1*(1 - (2*W)^m_net1)); 
funY=@(Z,W_net2,m_net2) (2*(1 - 2*Z))/((1 - 2*Z)*(1 + W_net2) + Z*W_net2*(1 - (2*Z)^m_net2));     
funW=@(X,Y,n_W,n_L) (1 - ((1 - X)^(n_W - 1))*(1 - Y)^n_L);                                            
funZ=@(X,Y,n_W,n_L)(1 - ((1 - X)^n_W)*(1 - Y)^(n_L - 1));                                             
funS=@(X,Y,n_W,n_L,beta)(Y*((1 - Y)^(n_L - 1))*(1 - X)^(n_W)-beta);
fun1=@(X,Y,W,Z,W_net2,m_net2,n_L,W_net1,m_net1,n_W,beta) ([funX(W,W_net1,m_net1); funY(Z,W_net2,m_net2) ; funW(X,Y,n_W,n_L); funZ(X,Y,n_W,n_L) ; funS(X,Y,n_W,n_L,beta)]); 

n_W = 3;
n_L = 4;
beta = 0.45; 
m_net2 = 4;
W_net1= 16;
m_net1= 6;
fun2=@(P) (P-fun1(P(1),P(2),P(3),P(4),P(5),m_net2,n_L,W_net1,m_net1,n_W,beta));
InitialGuess=[0;0;0;0;0];
fsolve(fun2,InitialGuess)

and run the code, I get

ans =

0.0059

1.1505

1.0010

1.0604

0.0354

When I use this values, they do not satisfy the given equations. For example with the given values of P(3)=W=1.0010, W_net1= 16,m_net1= 6 and by substituting them in the equation of X (funX), the value that I got is different from P(1) given by fsolve function. Do I do somthing wrongly?

3- Is it possible to add constraints to these equations? I mean is there any way that MATLAB knows 0<=P(1)-P(4) <=1 (they are probability).

Thank you so much in advance.

[Susan]

Thanks for the note.

The P(1)-P(4) and beta are probabilities and their values should be between zero and 1. Moreover, W_net2 \in {16,32,64,128,256,512,1024}.

Then according the tests, there is no solutions....

I am pretty sure that the equations are correct. So, how can I deal with this situation?

Thanks again

[Walter Roberson]

I took the same approach that you did for fun2: I made W_net2 the fifth element of P.

The easiest way to see the 5 equations is to use

P = sym('p', [5 1]);
fun2(P)

It is not possible to add constraints to fsolve(). It is, however, possible to add constraints to vpasolve():

fun2s = fun2(P);
sol = vpasolve(fun2s, [0 1; 0 1; 0 1; 0 1; -10 10]);  %adjust as appriopriate for W_net2

[Walter Roberson]

I see your note about W_net2 that restricts it to some powers of 2. I thought W_net2 was intended to be output?? It is very very unlikely that probabilities would lead to exact powers of 2 for W_net2: those are more like what you would use for input, perhaps solving the system once for each discrete W_net2 . (Which would be a problem because you would then have 5 equations in four unknowns.)

[Susan]

Thanks again for your reply. I have learnd a lot from you. Thanks!

You are rigth. W_net2 is an output. But the value I expect it takes is one of {16,32,64,128,256,512,1024} values or at least sth in this range.

BTW, do you know why I ended up with this situation that I explain above:

"When I use this values, they do not satisfy the given equations. For example with the given values of P(3)=W=1.0010, W_net1= 16,m_net1= 6 and by substituting them in the equation of X (funX), the value that I got is different from P(1) given by fsolve function. Do I do somthing wrongly? "

Thank you so much

[Walter Roberson]

fsolve() has a termination tolerance. In some cases, especially cases with exp() or with division, using a tighter tolerance can make a big difference after back substitution.

[Susan]

I see. Thank you very much for your response.

[Walter Roberson]

Your equations have singularities when P(3) = 1/2 or P(4) = 1/2. That is an odd thing to have when you are calculating with probabilities.

[Walter Roberson]

When I fsolve() with no options, it runs out of function evaluations. When I pass in options that permit more evaluations, it stops saying last step was ineffective but result is not near 0 -- in other words it is getting caught in a local minima. fsolve() is not a global root solver: it depends on the initial value to figure out where to proceed.

[Walter Roberson]

vpasolve() restricted to [0 1; 0 1; 0 1; 0 1; 0 2000] takes a while, but eventually returns empty, indicating that there are no solutions it could find in that range.

[Susan]

Thanks for the note.

I tried vpasolve() and it is still running. As you said, it takes a while. Good to hear that you have already tried that and it returend empty.

You brought up a very good point. The equations have singularities when P(3) = 1/2 or P(4) = 1/2. However, since both nominator and denominator have the value of 0 at the same time lim W and Z when p->1/2 is available. I don't know how matlab deals with these situations. Is it possible that this makes the problem?

Again, thank you SO much.

[Walter Roberson]

When working numerically you would get 0/0 which would be NaN.

The tests I am doing effectively recover from NaN, so I know that is not the reason we cannot find a root.

[Susan]

I see. Thanks for letting me know.

I am trying another approach to see whether it works out or not.

This is what I am doing now

1) knowing beta, from S ===> relationship between X and Y can be found. (using solve function)

2) from X and W ===> X and Y can be found numerically (perhaps fsolve. I am at this stage now)

3) based on Y, W_net2 can be derived

I don't know if it's going to give me any solution or not though.

Thanks again VERY much for all your help. I appreciate that

[Susan]

A silly question. If I have the following equations, have I defined the above functions correctly?

X = 2*(1 - 2*W))/((1 - 2*W)*(1 + W_net1) + W*W_net1*(1 - (2*W)^m_net1); 
Y = 2*(1 - 2*Z))/((1 - 2*Z)*(1 + W_net2) + Z*W_net2*(1 - (2*Z)^m_net2);     
W = 1 - ((1 - X)^(n_W - 1))*(1 - Y)^n_L;                                            
Z = 1 - ((1 - X)^n_W)*(1 - Y)^(n_L - 1);                                             
S = Y*((1 - Y)^(n_L - 1))*(1 - X)^(n_W);

that we know S is given and it is equal to beta.

By modifying the functions as follows (Please put me rigth if I am wrong) , P(1)-P(4) are n the desired range. However, W_net2 still gets negative value.

funX=@(X, W,   W_net1,   m_net1) (2*(1 - 2*W))/((1 - 2*W)*(1 + W_net1) + W*W_net1*(1 - (2*W)^m_net1) - X); 
funY=@(Y, Z,   W_net2,   m_net2) (2*(1 - 2*Z))/((1 - 2*Z)*(1 + W_net2) + Z*W_net2*(1 - (2*Z)^m_net2) - Y);     
funW=@(W, X,Y,n_W,n_L) (1 - ((1 - X)^(n_W - 1))*(1 - Y)^n_L - W);                                           
funZ=@(Z, X,Y,n_W,n_L)(1 - ((1 - X)^n_W)*(1 - Y)^(n_L - 1) - Z);                                             
funS=@(X,Y,n_W,n_L,beta)(Y*((1 - Y)^(n_L - 1))*(1 - X)^(n_W))-beta;

Thanks in advance

[Walter Roberson]

Your equations have bad brackets on X and Y. You have what appears to be an extra ) just before the / and you have what appears to be a missing ) at the end of them.

[Susan]

Got you. Thank you so much

[Walter Roberson]

With this newer set of equations I am able to get down to a long expression in W_net2. The expression is rather long, but turns out to be closed form, all "algebraic" numbers (it involves a cubic and a quartic.) Of the 12 combinations of expressions, 10 of them turn out to be complex valued in the range W_net2 = 0 to 1024. Two of the combinations are real valued -- but only over 0 to slightly less than 4. One of those two has a singularity in that range and so ends up not crossing 0. The other is a bit oddly shaped but does not cross zero in the range. So... there does not appear to be any real-valued solutions.

[Susan]

Thanks for the note.

Nice work! I appreciate your time. Just a quick question, is funS founction written correctly? shouldn't be "funS=@(X,Y,n_W,n_L,beta)(Y*((1 - Y)^(n_L - 1))*(1 - X)^(n_W) -beta);"? I mean beta is outside the () in the above equation.

Still I get p2>1 and W_net2<0 even with new equations.

[Susan]

Walter, Could you please kindly take a look at my code and tell me why I am not getting the results that you get? Thank you so much in advance.

clear;
clc;

funX=@(X, W, W_net1, m_net1) (   (2*(1 - 2*W)/((1 - 2*W)*(1 + W_net1) + W*W_net1*(1 - (2*W)^m_net1))) - X   );       
funY=@(Y, Z, W_net2 ,m_net2) (     (2*(1 - 2*Z)/((1 - 2*Z)*(1 + W_net2) + Z*W_net2*(1 - (2*Z)^m_net2))) - Y      );       
funW=@(W, X, Y, n_W, n_L) (    1 - ((1 - X)^(n_W - 1))*((1 - Y)^n_L) - W   );                                                 
funZ=@(Z, X, Y, n_W, n_L) (    1 - ((1 - X)^n_W)*((1 - Y)^(n_L - 1)) - Z   );                                                 
funS=@(X, Y, n_W, n_L, Ps) (  Y*((1 - Y)^(n_L - 1))*((1 - X)^n_W) - Ps  );                                         
fun1=@(X, Y, W, Z, W_net2, m_net2, n_L, W_net1, m_net1, n_W, Ps) ([funX(X, W, W_net1, m_net1); funY(Y, Z, W_net2, m_net2) ; funW(W, X, Y, n_W, n_L); funZ(Z, X, Y, n_W, n_L) ; funS(X, Y, n_W, n_L, Ps)]); 
n_W = 3;
n_L = 4;
Ps = 0.45; 
m_net2 = 4;
W_net1 = 16;
m_net1 = 6;
fun2=@(P) (P-fun1(P(1), P(2), P(3), P(4), P(5), m_net2, n_L, W_net1, m_net1, n_W, Ps));
InitialGuess=[0;0;0;0;16];
fsolve(fun2,InitialGuess)

P = sym('p', [5 1]);
fun2(P)

Answer 1

[John D'Errico]

syms a x b y
EQ = y == a*x + b;
isolate(EQ,x)
ans =
x == -(b - y)/a

[Susan]

Thank you so much John for your reply. I have got the following error:

Undefined function 'isolate' for input arguments of type 'sym'.

Any idea?

If you have any idea to how to deal with the second part of my question, I would love to hear. Thanks in advance

[John D'Errico]

Are you saying that you tried exactly the code I wrote, and it complained about isolate?

What version of MATLAB do you have? isolate was introduced before R2006a, so I think it unlikely that you do not have isolate. And if it recognizes that sym is a valid class, then you should have the symbolic toolbox. Do you?

If you do not have the symbolic TB, then you cannot have MATLAB do a symbolic computation as you asked. Of course, you have written this all using the fsolve as a target. So I need to know what you have. Are you really looking for a symbolic solution as you initially say, or are you looking for a numerial solution, or do you care as long as you find a solution?

I cannot easily go further to help you unless you can at least verify what you have, so I might know what to suggest.

[Susan]

Yes,

I use R2014a and don't know why I get this error.

I tried the following code and get the same error

syms x a b c

eqn = a*x^2 + b*x + c == 0;

xSol = isolate(eqn, x)

Update: Interesting! I tried the code on R2016b and get the same error. What's wrong?

>> syms a x b y

EQ = y == a*x + b;

isolate(EQ,x)

Undefined function or variable 'isolate'.

[John D'Errico]

Strange. R2018b here. In that quadratic equation, I get this:

syms x a b c
eqn = a*x^2 + b*x + c == 0;
xSol = isolate(eqn, x)
xSol =
x == -(b + (b^2 - 4*a*c)^(1/2))/(2*a)

Yet, doc isolate tells me this:

"Introduced before R2006a"

So it clearly has been around sufficiently long, and it seems unlikely that isolate has changed functionality.

Is the symbolc toolbox working for you otherwise? For example, does solve work?

solve(EQ,x)
ans =
-(b - y)/a

or

solve(eqn,x)
ans =
 -(b + (b^2 - 4*a*c)^(1/2))/(2*a)
 -(b - (b^2 - 4*a*c)^(1/2))/(2*a)

Have you got some other m-file that you named isolate?

which isolate -all
/Applications/MATLAB_R2018b.app/toolbox/symbolic/symbolic/@sym/isolate.m  % sym method

[Walter Roberson]

Is it possible that you have Maple installed on your system?

[Susan]

@John, thanks for your reply. I appreciate your help. I don't have any other m-file that named isolate. And regarding your question on symbolc toolbax, yes it works fine.

@Walter, How much is it going to help me on solving this problem? If it is really needed, I can buy that. Thanks again.

[Walter Roberson]

Su san: Maple offers a symbolic toolbox interface for MATLAB that differs in the available commands. When we see people with some commands working abut not other commands, sometimes it turns out that they are using the Maple interface in MATLAB instead of MATLAB's own Symbolic Toolbox. (But that turns out not to be the case this time.)

John: When I "doc isolate" I am given the documentation for syms . It takes more work to find the documentation for isolate. It turns out that it was introduced in R2017a as a MATLAB interface, which is after Su san's R2014a release.

[Susan]

@ Walter, I see. Thanks for the clarification.

[Susan]

I would like to extend the above equation to the following ones. Any idea how I can implement them. Thanks in advance.

for (k = 1: K), (i = 1 : Nw), and (j = 1: Nl)
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))); 
W(i,k) = 1 - prod_{ii = 1, ii ~= i}^{Nw} (1 - X(ii, k) ) * prod_{j = 1}^{Nl} (1 - Y(j, k));
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))); 
Z(j,k) = 1 - prod_{jj = 1, jj ~= j}^{Nl} (1 - Y(jj, k) ) * prod_{i = 1}^{Nw} (1 - X(i, k));
S(j,k) = Y(j, k) * prod_{jj = 1, jj ~= j}^{Nl} (1 - Y(jj, k) ) * prod_{i = 1}^{Nw} (1 - X(i, k));

[Walter Roberson]

ivec = 1 : Nw;
jvec = 1 : Nl;
X = zeros(Nw, K);
Y = zeros(Nl, K);
W = zeros(Nw, K);
Z = zeros(Nl, K);
S = zeros(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

[Susan]

Thank you so much Walter! I appreciate your help.

Could you please kindly tell me how I can define fun1 anf fun2 in this case?

Now the parameters that I am looking for are X (Nw, K), Y(Nl, K), W(Nw, K), Z(Nl, K) and W_net2 (Nl, K). In this case S (Nl,K) is given and fixed.

Thank you in advance. I appreciate your time.

[Walter Roberson]

Could you please kindly tell me how I can define fun1 and fun2 in this case?

.... ummm, No?

Everything before used scalar X, Y, W, Z and what you would mean with arrays for those is not clear.

You could vectorize all of those previous functions, changing all of the * to .* and all of the ^ to .^ to make calculations on entire arrays at once. On the other hand your arrays are not all the same size, so some of those * should perhaps be left as * in order to invoke algebraic matrix multiplication ??

But then when you get to the fsolve stage, it is not clear what you would be solving. If you were to ask to fsolve an array of Nw x Nl x 4 values, then fsolve would be trying to find parameters that would work for all of the entries simultaneously, which probably is not what would be desired.

[Susan]

Thanks for your reply and my appologies for not being clear.

Accualty, everything now use scalar too, as before. The only diference is now these scalar depends on k, ivec, jvec. For example, the previus fun1 and fun 2 were only for one k and the X for all Nw (Y for all Nl) was equal. Now, the value of Y for each n_l (\in jvec) would be diferent from the value of Y for n_lprime (\in jvec). The same for X. Does it make sence?

Assuming all the X (Y) are equal for ivec = 1 : Nw (jvec = 1 : Nl;), I think if we put the functions in a for loop that goes from 1 to K, we are able to find X(:,1: K), Y(:, 1:K) ,.... Am I right?

The dificult part for me is writhing Nw+Nl functions, since fun1 and fun2 are a functions of this Nw+Nl functions.

Is it clear now? Any suggestions would be greatly appreciated.

Thanks again.

[Susan]

Hi Guys,

Any suggestions on how to apply fsolve to these functions?

Thank you so much in advance.

[Walter Roberson]

Yes you can use nested for loops.

You might perhaps also be able to use ndgrid() to construct the input values, and then use arrayfun() to do the fsolve()

[Susan]

Thanks for your quick reply.

All these functions are new to me. Let me to do a quick search on what you suggested and see if I underestand what you are saying.

Thanks again for all your help. I DO appreciate that.

[Susan]

Walter,

Would you please be so kind as to help me to implement your suggestion in MATLAB? I am totally lost.

Thank you very much

[Walter Roberson]

I will not get to this right away; I have something else I need to do tonight.

[Susan]

I underestand. Whenever you get time works. I really appreciate your time and help. Thanks again for everything.

[Susan]

Walter, could you please help me to implement your suggestion whenever you have time? Thanks in advance.

Answer 2

[Susan]

I have got some silly questions. Sorry in advance if they are super simple. Could anyone kindly help me to figure them out?

1) As for a variable that I am looking for, MATLAB gives me

W_net2 = 1665452806855272103430873159617883565696871742382633107303767067803714653521912534085508640394758232435502244395851440048832512/13020781553015561664280552453950494568153170161082649667490158510277680133267873702693313295722419712822602179610396291511542125

How is it possible that I simply get 0.1279 instead of this very long result? Why doesn't MATLAB do the division?

2) When I use solve(), there is a "z" in results. According to my underestanding, by using vpa() I should get rid of z, but I don't. Any idea?

Thanks in advance.

[Walter Roberson]

1)

vpa(W_net2)

The Symbolic Toolbox is primarily for working with infinitely precise solutions, closed form when possible. By default it converts all floating point numbers (which are imprecise) to reasonably close algebraic numbers or multiples of Pi, and then works with those in order to try to find precise results.

It is possible to get the Symbolic Toolbox to work with software floating point numbers, but there are surprising limitations. Including that A + B - A might not exactly equal B which is a fundamental floating point limitation. (But the Symbolic Toolbox manipulates formulas under the assumption that it can work algebraically with association and commutation and distribution, so the answers it comes up with when you use software floating point might not be the same as what you might arrive at if you calculate by hand...)

2)

z can show up in three different contexts:

• it can show up in root(), in the form root(polynomial in z, z) or root(polynomial in z, z, number) . These stand in for the set of values, z, such that the polynomial evaluated at z becomes 0 -- the roots of the polynomial. z is a placeholder variable in this situation, that does not necessarily correspond to any of your variables, but one of your variables involves an expression that involves roots of a polynomial. root() was added in R2015b. vpa() of root() will resolve down to specific numeric values only in the case where all of the coefficients of the polynomial are numeric or standard functions applied to a numeric value (e.g., z^5 - cos(pi/7)*z + exp(2))
• it can show up in RootOf(), in the form RootOf(expression in z, z) or RootOf(expression in z, z, number). This is similar to root() but the expression does not have to be a polynomial. These stand in for the set of values, z, such that the expression evalulated at z becomes 0 -- the roots of the expression. z is a placeholder variable in this situation, that does not necessarily correspond to any of your variables, but one of your variables involves an expression that involves roots of an expression. RootOf are more likely to show up in older releases of MATLAB, but can still show up in cases where the expression is not a polynomial. vpa() of RootOf() is only certain to resolve down to specific numeric values only in cases similar to what root() can handle. Because RootOf() might involve more complicated functions such as erfc() or ilaplace or integrals, even when there are no other unbound symbolic variables, MATLAB might not be able to find a software floating point approximation, so vpa() of a RootOf() sometimes returns a RootOf() rather than a software floating point number.
• it can show up in the results of solve() outside of a root() or RootOf() . When this happens, you probably need to add the 'ReturnConditions', true option to solve() and assign the output of solve to a single variable, after which you should look at the Parameters field of the solution. If Parameters is non-empty, it will be a list of variables that the Symbolic Toolbox introduced to make the solution easier to write. For example if x and y both involve the result of calculating an expression, such as x = 3*expression^2 + 7, and y = expression/19, then instead of putting the full value of the expression into x and y, the Symbolic Toolbox might introduce a temporary variable, such as z, and then say x = 3*z^2 + 7, y = z/19 . You then need to look in the Conditions property of the solution to see how the temporary variable is defined or constrained. Typically the Symbolic Toolbox will only define a temporary variable if the variable can have multiple values, but it also sometimes defines a temporary variable for complicated conditions, such as "z has to be such that arccos(sqrt(z-2))/Pi is a positive integer and z = sin(a^3+log(b))" . Sometimes it is possible, with some work, to show that the conditions cannot be met. When a Parameter shows up in the result of solve(), then vpa() will never get rid of the parameter: you have to find solutions for the Conditions and subs() those into the solutions.

[Susan]

Thank you SO much for your through reply. It helps a lot.

[Susan]

Hey guys,

I have got some questions for you and really appreciate your help. I use the following code

n_L = 3; n_W = 2;  W_net1 = 16; m_net1 = 6; m_net2 = 6;  Beta_L = 0.25;
syms Y X
eqn = Y*((1-Y)^(nL - 1))*((1 - X)^nW) - Beta_L  == 0;
solY = solve(eqn, Y, 'MaxDegree', 4); % Add  'MaxDegree', 4 to get rid of z

"solY" give me [3 * 1] which each row is a function of X and this is what I am looking for. Then I use,

funX=@(X, W, W_net1, m_net1) (   (2*(1 - 2*W)/((1 - 2*W)*(1 + W_net1) + W*W_net1*(1 - (2*W)^m_net1))) - X   );       
funY=@(Y, X) (   solY(1,1) - Y      );       
funW=@(W, X, Y, n_W, n_L) (    1 - ((1 - X)^(n_W - 1))*((1 - Y)^n_L) - W   );                                                 
                                       
fun1=@(X, Y, W, W_net1, n_L, m_net1, n_W) ([funX(X, W, W_net1, m_net1); funY(Y, X) ; funW(W, X, Y, n_W, n_L)]); 

fun2=@(P) (P-fun1(P(1), P(2), P(3), W_net1, n_L, m_net1, n_W));
InitialGuess=[0;0;0];
fsolve(fun2,InitialGuess)

then I get the following error

Error using fsolve (line 269)

FSOLVE requires all values returned by functions to be of data type double.

Error in test (line 41)

x = fsolve(fun2,InitialGuess);

Anyone knows why? However, if I put the value of solY(1,1) on the funY function, I don't get any errors and the program runs perfectly.

Thanks in advance

[Walter Roberson]

Any relationship between n_L and nL ? Any relationship between n_W and nW ?

[Walter Roberson]

funY = @(Y, X) double( subs(solY(1,1), sym('X'), X) - Y);

[Susan]

Thank you so much Walter!

And sorry for the mess. Yes, n_L = nL and n_W = nW

Answer 3

[Susan]

Hey MATLAB experts,

Can anyone kindly help me to write down the functions in 'solve function' for below equations? Thanks in advance.

ivec = 1 : Nw;
jvec = 1 : Nl;
X = zeros(Nw, K);
Y = zeros(Nl, K);
W = zeros(Nw, K);
Z = zeros(Nl, K);
S = zeros(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

[Walter Roberson]

You have different sizes of arrays, so I do not know what is intended to correspond to what.

Do not try to be fancy: just use nested loops to find each combination of scalar parameters and run those through your solve() code and accumulate the output in an array of appropriate size. Get it working first before you worry about vectorzing it.

[Susan]

Thanks for your reply.

I wanted to use ndgrid() to construct the inputs and then use the arrayfun() as you suggested but unfortunately all my tries wasn't succesful.

Sure! I will use nested loop and try to get it working first. Thanks again for your reply.

[Susan]

Hi Walter,

I have tried to implement this in Matlab, but still get some errors. Would you please be so kind as to take a look at my code and tell me what kind of modification I need to do?Thanks in advance!

%%
clc
clear all
K = 5;
Nw = 2;
Nl = 3;
W_net1 = 16*ones(Nw, K);
m_net1 = 6*ones(Nw, K);
W_net2 = 16*ones(Nl, K);
m_net2 = 6*ones(Nl, K);
    
ivec = 1 : Nw;
jvec = 1 : Nl;
X = zeros(Nw, K);
Y = zeros(Nl, K);
W = zeros(Nw, K);
Z = zeros(Nl, K);
S = zeros(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     
                                                 
results = zeros(4, K);                                              
for k = 1 : K
    for i = ivec
        for j = jvec
            funX = @(X, W, W_net1, m_net1) (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));
            funY = @(Y, Z, W_net2 ,m_net2) (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));
            funW = @(W, X, Y, Nw, Nl) (1 - prod( 1 - X(setdiff(1 : Nw, i), k) ) * prod( 1 - Y(1 : Nl, k) ) - W(i,k));
            funZ = @(Z, X, Y, Nw, Nl) (1 -  prod(1 - Y(setdiff(1 : Nl, j), k)) * prod(1 - X(1 : Nw, k)) - Z(j,k));
            fun1=@(X, Y, W, Z, W_net1, m_net1, W_net2, m_net2, Nw, Nl) ([funX(X, W, W_net1, m_net1);...
                funY(Y, Z, W_net2 ,m_net2) ; funW(W, X, Y, Nw, Nl); funZ(Z, X, Y, Nw, Nl)]);
            
            fun2=@(P) (P-fun1(P(1), P(2), P(3), P(4), W_net1, m_net1, W_net2, m_net2, Nw, Nl));
            InitialGuess=[0;0;0;0];
            fsolve(fun2,InitialGuess)
            results(:, k) = [P(1); P(2); P(3); P(4)];
        end
    end
end

[Walter Roberson]

fsolve is going to invoke fun2 with a vector of length 4, received as variable P. The scalar P(1) is passed to fun1, where it becomes known as X. fun1 passes this scalar X to a number of calls, including as the second parameter to funW, where it is known as X (still). funW tries to index that X at a vector of locations.

[Walter Roberson]

Note that I am still lost as of https://www.mathworks.com/matlabcentral/answers/455589-solving-system-of-n-equations#comment_694341 . I have never recovered from that.

If you want to extend some of the variables to be arrays where they used to be scalars, then take the previous code and put it into a function that accepts one scalar at a time for those values, and call that function in nested loops to generate all combinations of input values that are relevant.

However, if you are extending to arrays with values corresponding to your previous X and Y positions, then I am lost about what the new equations are.

[Susan]

Thanks so much Walter for your comments.

I want to extend some of the variables to be arrays where they used to be scalars. As you suggested, took the previous code and put it into a function that accepts one scalar at a time for those values, i.e.,

function[P] = calculate_p(n_W,n_L,Wmin_LAA,m_LAA,Wmin_WiFi,m_WiFi)
funX=@(X, W, Wmin_WiFi, m_WiFi) ((2*(1 - 2*W)/((1 - 2*W)*(1 + Wmin_WiFi) + W*Wmin_WiFi*(1 - (2*W)^m_WiFi))) - X );        
funY=@(Y, Z, Wmin_LAA ,m_LAA) ( (2*(1 - 2*Z)/((1 - 2*Z)*(1 + Wmin_LAA) + Z*Wmin_LAA*(1 - (2*Z)^m_LAA))) - Y  );          
funW=@(W, X, Y, n_W, n_L) (    1 - ((1 - X)^(n_W - 1))*((1 - Y)^n_L) - W   );                                           
funZ=@(Z, X, Y, n_W, n_L) (    1 - ((1 - X)^n_W)*((1 - Y)^(n_L - 1)) - Z   );                                              
fun1=@(X, Y, W, Z, Wmin_LAA, m_LAA, n_L, Wmin_WiFi, m_WiFi, n_W) ([funX(X, W, Wmin_WiFi, m_WiFi); funY(Y, Z, Wmin_LAA, m_LAA) ; funW(W, X, Y, n_W, n_L); funZ(Z, X, Y, n_W, n_L) ]); 

    

       
fun2=@(P) (P-fun1(P(1), P(2), P(3), P(4), Wmin_LAA, m_LAA, n_L, Wmin_WiFi, m_WiFi, n_W));
InitialGuess=[0;0;0;0];
fsolve(fun2,InitialGuess) 
end

Now, as you mentioned, I have to call that function in nested loops to generate all combinations of input values that are relevant. for all i, j, k.

I am strugling in this point. where should I call the above function in the following one? SHould I prelocate X,Y,W, Z with zeros or it's going to give me some wrong results?

ivec = 1 : Nw;
jvec = 1 : Nl;
%X = zeros(Nw, K);
%Y = zeros(Nl, K);
%W = zeros(Nw, K);
%Z = zeros(Nl, K);
%S = zeros(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 

What I want to do is for a given K, Nw, Nl, solve the following equations simualtinously and figure out what is the value of X(i, k), W(i,k), Y(j, k) , Z(j,k) for each i, j, k.

For example if K = 1, Nw = 2, Nl = 3, I need to know what are the values of X(1,1), X(2, 1), Y(1,1), Y(2,1), Y(3,1), W(1,1), W(2, 1), Z(1,1), Z(2,1), and Z(3,1).

Any help would be greatly appreciated. Thanks again

[Susan]

Walter,

Could you please kindly help me to solve this issue?

Many thanks in advance

[Walter Roberson]

Sorry, I do no know when I will be well enough to address this.

[Susan]

No worries. Get well soon :)

Answer 4

[Susan]

Hi Walter,

In one of the comments above you mentioned that "The tests I am doing effectively recover from NaN, so I know that is not the reason we cannot find a root." Did you do any specific things that the test recover from Nan?

The reason I am asking is that I am using fmincon() to solve two optimization problems (the objective function is the same but I optimize the objfun w.r.t. two different variables), and regardless of what I am selecting the initial values, I ended up with this error

"Error using sqpInterface

Objective function is undefined at initial point. Fmincon cannot continue."

Thanks in advance.

[Walter Roberson]

I am not using fmincon(); I am using my own algorithm that includes a bunch of fminsearch() calls that I have wrapped around with some range restraints, convincing it to do things it was not designed to do. My wrappers start from an initial point and keep searching until the location is out of range or the function call returns non-finite, and when they find an undesired location, they back off to the last valid place they looked at and return that. (They do not try to "bounce off" of the invalid region though.)

When you use fmincon with sqp, your initial point will be manipulated if necessary so that it is at upper or lower bound instead of being outside of range. However, the initial point will not be manipulated to be within linear inequalities (not for sqp -- the situation is slightly different for trust-region.) It is best to start with a point that is within the bounds constraints and the linear inequalities, and preferably the linear equalities as well. It is common to not know an initial point that is within the nonlinear constraints.

When you get told the objective function is undefined at the initial point, it means you are getting nan or one of the infinities. Your objective function will be called with the input initial point (modified for upper/lower constriants if need be) so you can test that directly without worrying about the other constraints.

[Susan]

Thanks for your reply. It totaly makes sense.

I start with a point that is within the bounds constraints and don't have any linear inequalities. However, I have a complex nonlinear constraint.

Do you know how I can replace this NaNs with a value in a cell array? I have the following and would like to replace the r with a value when it is Nan.

r = cell(sum(nL), numel(nL), numel(nW), max(nW(:)));
for k = 1 : numel(nW)
    for m = 1 : nW(k)
        for j = 1 : numel(nL)
            for i = 1 : nL(j)
                r{i + Nup*(j - 1), j, k, m} = .....
              
            end
        end
    end
end
 %%
 r(cellfun(@(cell) any(isnan(cell(:))),r))={0}

(I have added a line after %%, is it a correct way to do so?)

Thank you so much in advance

[Susan]

It just came to my attention that I should use the tools explained on this page to make my question more readable. Just wanted to THANK YOU for being SO nice to edit all of my questions without even mentioning that! Thanks again Walter!

[Walter Roberson]

That would work, but I would recommend against using the variable named cell due to its use as the name of the constructor functions for cell arrays.

[Susan]

Thank you! I will use another name for this variable then! Thanks again