Fsolve
in [Matlab]
|
From: rabbahs on 20 Mar 2010 04:26 Hi all, i am trying to solve a system of linear/nonlinear equations using fsolve command in matlab. But i am getting some undesirable results. First of all, kindly look at the code. function F=myfun(X,c1) a=X(1); b=X(2); c=X(3); d=X(4); e=X(5); f=X(6); g=X(7); p=X(8); F(1)=a+b+d+g-c1(4); F(2)=c1(6)+2*c1(7)+c1(8)-2*a-b-f; F(3)=c1(5)+2*c1(8)-2*c-4*d-2*f; F(4)=3.76*c1(7)-e; F(5)=c1(1)*a*p-(b*b); F(6)=c1(2)*f*p-(b*c); F(7)=c1(3)*c*c-(d*p); F(8)=a+b+c+d+e+f-p; as u can see from above function script, i have 8 unknowns (i.e a,b,c,d,e,f,g & p). c1(1) to c1(8) are constant (having value greater then 0) the fsolve command which i am using to solve the above function is, fsolve(@(X) myfun_1(X,[q r s m x y z k]),[1 1 1 1 1 1 1 1],options); where, myfun_1 is function file name q r s m x y z k are constant c1(1) to c1(8) (all greater than zero) options = optimset('MaxFunEvals',3000,'MaxIter',20000); [11111111] are initial guess the problems i am facing with this code are, 1) most of the time the solution is not converged it gives me msgs like, change the initial guess ( which i changed thousands of time.) it ask me to increase the MaxFunEvals and iteration (it tried it to up to millions) 2) I am expecting values from this code between 0-5 for the variable (surely greater than zero) but it is giving me negative results. so is some thing wrong with my code ? is there any cure for it ? kindly let me know any possible solution for this problem, if any other information is required then kindly let me know. thanks and regard have a nice day bye.
From: Walter Roberson on 20 Mar 2010 13:38 rabbahs wrote: > Hi all, > > i am trying to solve a system of linear/nonlinear equations using fsolve > command in matlab. But i am getting some undesirable results. > > First of all, kindly look at the code. > > function F=myfun(X,c1) > > a=X(1); > b=X(2); > c=X(3); > d=X(4); > e=X(5); > f=X(6); > g=X(7); > p=X(8); > > F(1)=a+b+d+g-c1(4); > F(2)=c1(6)+2*c1(7)+c1(8)-2*a-b-f; > F(3)=c1(5)+2*c1(8)-2*c-4*d-2*f; > F(4)=3.76*c1(7)-e; > F(5)=c1(1)*a*p-(b*b); > F(6)=c1(2)*f*p-(b*c); > F(7)=c1(3)*c*c-(d*p); > F(8)=a+b+c+d+e+f-p; > > as u can see from above function script, i have 8 unknowns (i.e > a,b,c,d,e,f,g & p). > c1(1) to c1(8) are constant (having value greater then 0) Exact solution: Let P = RootOf( 10000*(4*c13 + 1)*(c11*c12^2*c13 + 4*c12^2*c13 - c11*c12 + c12^2 - c11)*_Z^4 + 800*(25*c11*c12^2*c13*c16 + 169*c11*c12^2*c13*c17 - 188*c11*c12*c13*c17 + 100*c12^2*c13*c16 + 1052*c12^2*c13*c17 - 25*c11*c12*c16 - 216*c11*c12*c17 + 25*c12^2*c16 + 263*c12^2*c17 - 25*c11*c16 - 169*c11*c17)*_Z^3 - 4*(625*c11*c12^2*c13*c15^2 + 2500*c11*c12^2*c13*c15*c16 + 16900*c11*c12^2*c13*c15*c17 + 2500*c11*c12^2*c13*c15*c18 + 5000*c11*c12^2*c13*c16*c18 + 33800*c11*c12^2*c13*c17*c18 + 2500*c11*c12^2*c13*c18^2 + 5000*c12^2*c13*c15^2 + 10000*c12^2*c13*c15*c16 + 105200*c12^2*c13*c15*c17 + 20000*c12^2*c13*c15*c18 + 20000*c12^2*c13*c16*c18 + 210400*c12^2*c13*c17*c18 + 20000*c12^2*c13*c18^2 - 625*c11*c12*c15^2 - 1250*c11*c12*c15*c16 - 13150*c11*c12*c15*c17 - 2500*c11*c12*c15*c18 + 2500*c11*c12*c16^2 + 43200*c11*c12*c16*c17 - 2500*c11*c12*c16*c18 + 177788*c11*c12*c17^2 - 26300*c11*c12*c17*c18 - 2500*c11*c12*c18^2 + 1250*c12^2*c15^2 + 2500*c12^2*c15*c16 + 26300*c12^2*c15*c17 + 5000*c12^2*c15*c18 - 2500*c12^2*c16^2 - 52600*c12^2*c16*c17 + 5000*c12^2*c16*c18 - 276676*c12^2*c17^2 + 52600*c12^2*c17*c18 + 5000*c12^2*c18^2 - 625*c11*c15^2 - 9400*c11*c15*c17 - 2500*c11*c15*c18 + 2500*c11*c16^2 + 33800*c11*c16*c17 + 78900*c11*c17^2 - 18800*c11*c17*c18 - 2500*c11*c18^2)*_Z^2 + 4*c12*(c15+2*c18)*(25*c15 + 50*c16 + 526*c17 + 50*c18)*(25*c11*c16 + 169*c11*c17 - 50*c12*c16 - 526*c12*c17)*_Z + 625*c12^2*c15^4 + 2500*c12^2*c15^3*c16 + 26300*c12^2*c15^3*c17 + 5000*c12^2*c15^3*c18 + 2500*c12^2*c15^2*c16^2 + 52600*c12^2*c15^2*c16*c17 + 15000*c12^2*c15^2*c16*c18 + 276676*c12^2*c15^2*c17^2 + 157800*c12^2*c15^2*c17*c18 + 15000*c12^2*c15^2*c18^2 + 10000*c12^2*c15*c16^2*c18 + 210400*c12^2*c15*c16*c17*c18 + 30000*c12^2*c15*c16*c18^2 + 1106704*c12^2*c15*c17^2*c18 + 315600*c12^2*c15*c17*c18^2 + 20000*c12^2*c15*c18^3 + 10000*c12^2*c16^2*c18^2 + 210400*c12^2*c16*c17*c18^2 + 20000*c12^2*c16*c18^3 + 1106704*c12^2*c17^2*c18^2 + 210400*c12^2*c17*c18^3 + 10000*c12^2*c18^4 ) The RootOf() operator means to find the values of _Z that make the expression 0. Notice that this is a quartic expression in _Z and thus has exact solutions -- but printing out the exact solutions would be rather long! Easier to substitute in your known constants and use Matlab's roots() operator to get the roots -- some of which may well be complex. And you might need the complex roots in order to proceed further. Once P are determined, Let Q = 50*(2*c12*c13 - 1)*P - 25*c15 - 50*c16 - 526*c17-50*c18 Let R = 2*(c12 + 2)*P - c12*c15 - 2*c12*c18 Then: a = (2500*(4*c13+1)*(c12^2*c13 - c12 - 1)*P^3 + (5000*c12^2*c13*c16 + 33800*c12^2*c13*c17 - 37600*c12*c13*c17 - 5000*c12*c16 - 43200*c12*c17 - 5000*c16 - 33800*c17)*P^2 + (-625*c12^2*c13*c15^2 - 2500*c12^2*c13*c15*c16 - 16900*c12^2*c13*c15*c17 - 2500*c12^2*c13*c15*c18 - 5000*c12^2*c13*c16*c18 - 33800*c12^2*c13*c17*c18 - 2500*c12^2*c13*c18^2 + 625*c12*c15^2 + 1250*c12*c15*c16 + 13150*c12*c15*c17 + 2500*c12*c15*c18 - 2500*c12*c16^2 - 43200*c12*c16*c17 + 2500*c12*c16*c18 - 177788*c12*c17^2 + 26300*c12*c17*c18 + 2500*c12*c18^2 + 625*c15^2 + 9400*c15*c17 + 2500*c15*c18 - 2500*c16^2 - 33800*c16*c17 - 78900*c17^2 + 18800*c17*c18 + 2500*c18^2)*P + 625*c12*c15^2*c16 + 4225*c12*c15^2*c17 + 1250*c12*c15*c16^2 + 21600*c12*c15*c16*c17 + 2500*c12*c15*c16*c18 + 88894*c12*c15*c17^2 + 16900*c12*c15*c17*c18 + 2500*c12*c16^2*c18 + 43200*c12*c16*c17*c18 + 2500*c12*c16*c18^2 + 177788*c12*c17^2*c18 + 16900*c12*c17*c18^2) / (25*Q*R) b = ( -100*c12*(4*c13 + 1)*P^2 - 4*c12*(263*c17 + 25*c16)*P + c12*(c15+2*c18)*(25*c15+50*c16 + 526*c17 + 50*c18) ) / (50*R) c = P d = -25*c13*P*R/Q e = 94/25*c17 f = 1/2* (100*(4*c13 + 1)*P^2 + 4*(25*c16 + 263*c17)*P - 25*c15^2 - 50*c15*c16 - 526*c15*c17 - 100*c15*c18 - 100*c16*c18 - 1052*c17*c18 - 100*c18^2) / ((100*c12*c13 - 50)*P - 25*c15 - 50*c16 - 526*c17 - 50*c18) g = ( 5000*(4*c12^2*c13^2 + 2*c12^2*c13 + 4*c12*c13 + 8*c13 + 1)*P^3 + 100*(100*c12^2*c13*c14 - 50*c12^2*c13*c15 + 376*c12^2*c13*c17 - 100*c12^2*c13*c18 + 200*c12*c13*c14 - 200*c12*c13*c15 - 200*c12*c13*c16 - 1352*c12*c13*c17 - 400*c12*c13*c18 - 50*c12*c14-25*c12*c15 - 188*c12*c17 - 50*c12*c18 - 100*c14+100*c16 + 676*c17)*P^2 - 2*(2500*c12^2*c13*c14*c15 + 5000*c12^2*c13*c14*c18 + 9400*c12^2*c13*c15*c17 + 18800*c12^2*c13*c17*c18 + 2500*c12*c14*c16 + 26300*c12*c14*c17 + 1250*c12*c15*c16 + 13150*c12*c15*c17 + 9400*c12*c16*c17 + 2500*c12*c16*c18 + 98888*c12*c17^2 + 26300*c12*c17*c18 + 2500*c14*c15 + 5000*c14*c16 + 52600*c14*c17 + 5000*c14*c18 + 625*c15^2 + 9400*c15*c17 + 2500*c15*c18 - 2500*c16^2 - 33800*c16*c17 - 78900*c17^2 + 18800*c17*c18 + 2500*c18^2)*P + 9400*c12*c15*c16*c17 + 71400*c12*c15*c17*c18 + 5000*c12*c15*c16*c18 + 18800*c12*c16*c17*c18 + 625*c12*c15^3 + 5000*c12*c18^3 + 98888*c12*c15*c17^2 + 17850*c12*c15^2*c17 + 3750*c12*c15^2*c18 + 7500*c12*c15*c18^2 + 1250*c12*c15^2*c16 + 5000*c12*c16*c18^2 + 197776*c12*c17^2*c18 + 71400*c12*c17*c18^2 + 26300*c12*c14*c15*c17 + 5000*c12*c14*c15*c18 + 5000*c12*c14*c16*c18 + 52600*c12*c14*c17*c18 + 1250*c12*c14*c15^2 + 5000*c12*c14*c18^2 + 2500*c12*c14*c15*c16 ) / (50*Q*R) p = -1/25*P*Q / ( (2*c12 + 4)*P - c12*c15 - 2*c12*c18 ) If you examine b, d, and p, you can see that it wouldn't take much for the values to come out negative. Substituting in the random constants [5/2, 9/5, 33/10, 31/10, 7/5, 17/10, 17/5, 23/10], I get the results: [a = 1.267103339, b = 9.012979781, c = 1.806899943, d = .4201432576, e = 12.78400000, f = .3528135419, g = -7.600226377, p = 25.64393986] [a = 18.20425353, b = -21.83226313, c = 2.310915533, d = 1.682664193, e = 12.78400000, f = -2.676243918, g = 5.045345414, p = 10.47332620] [a = 14.20536180, b = -18.69989687, c = -2.074911341, d = 1.442869037, e = 12.78400000, f = 2.189173266, g = 6.151666031, p = 9.846595894] [a = 1.749657476, b = 9.856441847, c = -5.904572379, d = 5.180164589, e = 12.78400000, f = -1.455756799, g = -13.68626391, p = 22.20993473] Notice that no matter which root of the quartic was chosen, at least one of the variables came out negative. I do not know at the moment if that is *necessarily* the case -- my computation engine has pretty much locked itself up processing the testing. I can say this: g is a function of c1(8) to the 12th power, involving quantities ranging from 10^(-10) to 10^54. With that range of values, order of operations would become crucial as round-off error is likely to *very* high. Scanning through the expression, if you were to calculate it naively in double precision, I can see that you would be lucky to get a single bit of the solution right. This high reliance upon the exact value of c18 is hidden in the way I rewrote the solution, above: P, the RootOf, involves up to c18^4, and g depends upon P^3 and thus upon (c18^4)^3 = c18^12 .
From: rabbahs on 20 Mar 2010 17:04 Mr. Walter, I really appreciate your effort to put your time to solve my problem. thanks what i get from you reply is that you find out the exact solution for the variable using Rootof( command in Matlab (i hope this is the case). what are these random constant ? c11 c12 ... as you find out that one of the value among variable always comes negative (does it mean that there is some problem in the equations it self, thats y it is giving negative values ?) how can i write the nonlinear equations in AX=B form. kindly refer me to some tutorial for fsolve (other then Matlab tutorial) again thanks have a nice day bye
From: Walter Roberson on 20 Mar 2010 19:20 rabbahs wrote: > how can i write the nonlinear equations in AX=B form. You can't. Some of your variables appear both squared and to the first power, so you cannot handle the situation by variable substitution. > what i get from you reply is that you find out the exact solution for > the variable using Rootof( command in Matlab (i hope this is the case). The Matlab equivalent of Maple's RootOf is roots() > what are these random constant ? c11 c12 ... c11 is what you wrote as c1(1) c12 is what you wrote as c1(2) and so on. I needed to rename during my processing because otherwise my symbolic engine would have thought that I was indicating function calls. > as you find out that one of the value among variable always comes > negative (does it mean that there is some problem in the equations it > self, thats y it is giving negative values ?) I discovered a moment ago that I had a typo when I input your equations originally, so the solution I gave is not valid!! I'm working towards a more compact solution now... or at least a more correct solution.
From: rabbahs on 21 Mar 2010 00:24
Mr. Walter, I really appreciate your effort to put your time to solve my problem. thanks what i get from you reply is that you find out the exact solution for the variable using Rootof( command in Matlab (i hope this is the case). what are these random constant ? c11 c12 ... as you find out that one of the value among variable always comes negative (does it mean that there is some problem in the equations it self, thats y it is giving negative values ?) how can i write the nonlinear equations in AX=B form. kindly refer me to some tutorial for fsolve (other then Matlab tutorial) again thanks have a nice day bye |