Solving system of n equations

Question

Susan 2019-4-10

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/455589-solving-system-of-n-equations

评论： Susan 2019-5-10

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

27 个评论
显示 25更早的评论隐藏 25更早的评论

Walter Roberson 2019-4-10

在 MATLAB Online 中打开

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.

Walter Roberson 2019-4-11

编辑：Walter Roberson 2019-4-11

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

Susan 2019-4-11

在 MATLAB Online 中打开

@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 2019-4-11

编辑：Walter Roberson 2019-4-11

在 MATLAB Online 中打开

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

Susan 2019-4-12

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 2019-4-12

在 MATLAB Online 中打开

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 2019-4-10

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/455589-solving-system-of-n-equations#answer_370070

在 MATLAB Online 中打开

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

20 个评论
显示 18更早的评论隐藏 18更早的评论

Susan 2019-4-16

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

Susan 2019-4-18

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

请先登录，再进行评论。

Answer 2

Susan 2019-4-12

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/455589-solving-system-of-n-equations#answer_370408

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.

6 个评论
显示 4更早的评论隐藏 4更早的评论

Walter Roberson 2019-4-12

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 2019-4-17

在 MATLAB Online 中打开

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 2019-4-17

在 MATLAB Online 中打开

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

Susan 2019-4-17

Thank you so much Walter!

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

请先登录，再进行评论。

Answer 3

Susan 2019-4-22

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/455589-solving-system-of-n-equations#answer_371697

在 MATLAB Online 中打开

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

9 个评论
显示 7更早的评论隐藏 7更早的评论

Susan 2019-5-8

编辑：Susan 2019-5-8

在 MATLAB Online 中打开

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

Susan 2019-5-9

编辑：Susan 2019-5-10

在 MATLAB Online 中打开

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

Walter Roberson 2019-5-10

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

Susan 2019-5-10

No worries. Get well soon :)

请先登录，再进行评论。

Answer 4

Susan 2019-4-26

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/455589-solving-system-of-n-equations#answer_372491

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.

5 个评论
显示 3更早的评论隐藏 3更早的评论

Walter Roberson 2019-4-26

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.

Walter Roberson 2019-4-26

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 2019-4-26

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

请先登录，再进行评论。

Solving system of n equations

27 个评论
显示 25更早的评论隐藏 25更早的评论

回答（4 个）

20 个评论
显示 18更早的评论隐藏 18更早的评论

6 个评论
显示 4更早的评论隐藏 4更早的评论

9 个评论
显示 7更早的评论隐藏 7更早的评论

5 个评论
显示 3更早的评论隐藏 3更早的评论

另请参阅

类别

标签

Community Treasure Hunt

Solving system of n equations

27 个评论 显示 25更早的评论隐藏 25更早的评论

回答（4 个）

20 个评论 显示 18更早的评论隐藏 18更早的评论

6 个评论 显示 4更早的评论隐藏 4更早的评论

9 个评论 显示 7更早的评论隐藏 7更早的评论

5 个评论 显示 3更早的评论隐藏 3更早的评论

27 个评论
显示 25更早的评论隐藏 25更早的评论

20 个评论
显示 18更早的评论隐藏 18更早的评论

6 个评论
显示 4更早的评论隐藏 4更早的评论

9 个评论
显示 7更早的评论隐藏 7更早的评论

5 个评论
显示 3更早的评论隐藏 3更早的评论