Solving System of Highly Non-Linear Equations
Ashley
syms p1 p2 p3 p4 p5 p6 epp2
eq1 = '((x1(i-1)/x1(i))^(p1-1))*dx1(i) + epp2*(x1(i-1)^p1-x1(i)^p1)*(x1(i-1)^(p1-1))*p1=dx1(i-1)';
eq2 = '((x2(i-1)/x2(i))^(p2-1))*dx2(i) + epp2*(x2(i-1)^p2-x2(i)^p2)*(x2(i-1)^(p2-1))*p2=dx2(i-1)';
eq3 = '((x3(i-1)/x3(i))^(p3-1))*dx3(i) + epp2*(x3(i-1)^p3-x3(i)^p3)*(x3(i-1)^(p3-1))*p3=dx3(i-1)';
eq4 = '((x4(i-1)/x4(i))^(p4-1))*dx4(i) + epp2*(x4(i-1)^p4-x4(i)^p4)*(x4(i-1)^(p4-1))*p4=dx4(i-1)';
eq5 = '((x5(i-1)/x5(i))^(p5-1))*dx5(i) + epp2*(x5(i-1)^p5-x5(i)^p5)*(x5(i-1)^(p5-1))*p5=dx5(i-1)';
eq6 = '((x6(i-1)/x6(i))^(p6-1))*dx6(i) + epp2*(x6(i-1)^p6-x6(i)^p6)*(x6(i-1)^(p6-1))*p6=dx6(i-1)';
eq7_1 = 'gx(i)';
eq7_2=' + (dx1(i)*((x1(i)^(1-p1))/p1)*(x1(i-1)^p1-x1(i)^p1)) + (dx2(i)*((x2(i)^(1-p2))/p2)*(x2(i-1)^p2-x2(i)^p2)) + (dx3(i)*((x3(i)^(1-p3))/p3)*(x3(i-1)^p3-x3(i)^p3)) + (dx4(i)*((x4(i)^(1-p4))/p4)*(x4(i-1)^p4-x4(i)^p4)) + (dx5(i)*((x5(i)^(1-p5))/p5)*(x5(i-1)^p5-x5(i)^p5)) + (dx6(i)*((x6(i)^(1-p6))/p6)*(x6(i-1)^p6-x6(i)^p6))';
eq7_3=' + 0.5*epp2*( ((x1(i-1)^p1-x1(i)^p1)^2) + ((x2(i-1)^p2-x2(i)^p2)^2) + ((x3(i-1)^p3-x3(i)^p3)^2) + ((x4(i-1)^p4-x4(i)^p4)^2) + ((x5(i-1)^p5-x5(i)^p5)^2) + ((x6(i-1)^p6-x6(i)^p6)^2) )= gx(i-1)';
eq7=sprintf('%s%s%s',eq7_1,eq7_2,eq7_3);
I am trying to use:
[p1,p2,p3,p4,p5,p6,epp2]=solve(eq1,eq2,eq3,eq4,eq5,eq6,eq7)
But is keeps telling me:
Warning: Explicit solution could not be found.
> In solve at 98
Is there a better way to approach this problem?
Thank you for your time,
Ashley
Torsten Hennig
I think you only have a chance to solve your system
numerically. Try fsolve.
Best wishes
Torsten.
Ashley
fsolve does not seem to be able to solve this problem. fsolve will not even iter and outputs the point x0 for the variables that I am trying to solve for.
Does anyone know other anyway to do this?
Torsten Hennig
fsolve does not accept symbolic variables.
Maybe that's the reason why it fails ?
Best wishes
Torsten.
Ashley
At the moment I cannot get fsolve to even run. It keeps saying:
??? Undefined function or method 'myfun' for input arguments of type 'double'.
Error in ==> @(x)myfun(gX,x1,x2,dx1,dx2)
I have no idea why I am getting this error. Here is what I coded:
x1=[40;50];
x2=[29.6645;47.9194];
gX=[600;21.5041];
dx1=[50;47.9194];
dx2=[40;29.6645];
X0 = [1;1;1]; % Make a starting guess at the solution
options=optimset('Display','iter');
[X,fval]=fsolve(@(x)myfun(gX,x1,x2,dx1,dx2),X0,options);
p1=x(1)
p2=x(2)
e=x(3)
In separate file: (file name myfun.m)
function F = myfun(gX,x1,x2,dx1,dx2)
F = [((x1(1)/x1(2))^(x(1)-1))*dx1(2) + x(3)*(x1(1)^x(1)-x1(2)^x(1))*(x1(1)^(x(1)-1))*x(1)+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0-dx1(1);...
((x2(1)/x2(2))^(x(2)-1))*dx2(2) + x(3)*(x2(1)^x(2)-x2(2)^x(2))*(x2(1)^(x(2)-1))*x(2)+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0-dx1(2);...
gX(2) + (dx1(2)*((x1(2)^(1-x(1)))/x(1))*(x1(1)^x(1)-x1(2)^x(1))) + (dx2(2)*((x2(2)^(1-x(2)))/x(2))*(x2(1)^x(2)-x2(2)^x(2))) + 0.5*x(3)*( ((x1(1)^x(1)-x1(2)^x(1))^2) + ((x2(1)^x(2)-x2(2)^x(2))^2) )-gX(1)];
The repeated ‘+0’ is so that the each row of the matrix has the same number of chars.
I realize looking at the code is a pain, but I would be very grateful for any help you can provide as I am really stuck.
Ashley
Thank you for your time,
Ashley
Miroslav Balda
> I have simplified my problem. I am now only trying to solve for three variables and I’m no longer using smys variables Here is what I have:
SNIP
Hi Ashley ,
I have tested your example by my function LMFnlsq from
www.mathworks.com/matlabcentral/fileexchange/17534
Your code has been changed into:
% Ashley.m 2009-11-11
global x1 x2 gX dx1 dx2
clc
x1=[40;50];
x2=[29.6645;47.9194];
gX=[600;21.5041];
dx1=[50;47.9194];
dx2=[40;29.6645];
X0 = [1,1,1]; % Make a starting guess at the solution
[x,ssq,cnt] = LMFnlsq('myfun',X0,'Display',-200,'maxiter',1605)
p1=x(1)
p2=x(2)
e=x(3)
The function defining residuals is:
function F = myfun(x)
global x1 x2 gX dx1 dx2
F = [((x1(1)/x1(2))^(x(1)-1))*dx1(2)+x(3)*(x1(1)^x(1)-x1(2)^x(1))*(x1(1)^(x(1)-1))*x(1)-dx1(1)
((x2(1)/x2(2))^(x(2)-1))*dx2(2)+x(3)*(x2(1)^x(2)-x2(2)^x(2))*(x2(1)^(x(2)-1))*x(2)-dx1(2)
gX(2)+(dx1(2)*((x1(2)^(1-x(1)))/x(1))*(x1(1)^x(1)-x1(2)^x(1)))+(dx2(2)*((x2(2)^(1-x(2)))/x(2))*(x2(1)^x(2)-x2(2)^x(2)))+0.5*x(3)*(((x1(1)^x(1)-x1(2)^x(1))^2)+((x2(1)^x(2)-x2(2)^x(2))^2))-gX(1)
];
The answer obtained by running it has been as follows:
******************************************************************
itr nfJ SUM(r^2) x dx
******************************************************************
0 1 1.9130e+006 1.0000e+000 0.0000e+000
1.0000e+000 0.0000e+000
1.0000e+000 0.0000e+000
200 286 2.0951e+004 1.8130e+000 -9.2236e-003
1.8655e+000 -1.0078e-002
3.5806e-003 3.1659e-004
400 576 1.6740e+004 2.3843e+000 -1.6110e-003
2.4954e+000 -1.9524e-003
1.9704e-005 3.1413e-007
600 884 1.4544e+004 2.7504e+000 -1.1610e-005
2.9067e+000 -1.6026e-005
7.0593e-007 1.8422e-009
800 1186 1.2638e+004 3.1171e+000 -2.8428e-005
3.3312e+000 -3.7778e-005
2.3648e-008 7.0455e-012
1000 1443 9.1712e+003 3.7514e+000 -2.2172e-004
4.3626e+000 -5.7972e-003
6.8743e-012 2.9561e-013
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 4.161889e-017.
> In LMFnlsq at 303
In Ashley at 12
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 4.161889e-017.
> In LMFnlsq at 303
In Ashley at 12
1155 1605 6.4117e+003 3.2964e+000 -5.8290e-004
5.6418e+000 -2.2968e-004
2.9771e-016 -5.1603e-019
x =
3.2964
5.6418
0.0000
ssq =
6.4117e+003
cnt =
-1605
p1 =
3.2964
p2 =
5.6418
e =
2.9771e-016
------------------------------------------------------
The option maxiter has been chosen equal 1605 in order that you see that your problem is ill-conditioned, even that the sum of squares dropped down from 1.9e6 to 6.4e3. That the solution is not final is also seen from the negative value of cnt, what indicates, that iterations have not been finished in maxiter iterations. It is all, what can I tell you at a moment.
Best regards,
Mira
Torsten Hennig
> using smys variables Here is what I have:
> At the moment I cannot get fsolve to even run. It
> keeps saying:
>
> ??? Undefined function or method 'myfun' for input
> arguments of type 'double'.
> Error in ==> @(x)myfun(gX,x1,x2,dx1,dx2)
>
> I have no idea why I am getting this error. Here is
> what I coded:
>
> x1=[40;50];
> x2=[29.6645;47.9194];
> gX=[600;21.5041];
> dx1=[50;47.9194];
> dx2=[40;29.6645];
>
> X0 = [1;1;1]; % Make a starting guess at the solution
> options=optimset('Display','iter');
> [X,fval]=fsolve(@(x)myfun(gX,x1,x2,dx1,dx2),X0,options
> );
[X,fval] = fsolve(@(x)myfun(x,gX,x1,x2,dx1,dx2),X0,options);
> p1=x(1)
> p2=x(2)
> e=x(3)
>
> In separate file: (file name myfun.m)
>
> function F = myfun(gX,x1,x2,dx1,dx2)
function F = myfun(x,gX,x1,x2,dx1,dx2)
> F = [((x1(1)/x1(2))^(x(1)-1))*dx1(2) +
> x(3)*(x1(1)^x(1)-x1(2)^x(1))*(x1(1)^(x(1)-1))*x(1)+0+0
> +0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0
> +0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0
> -dx1(1);...
> ((x2(1)/x2(2))^(x(2)-1))*dx2(2) +
> x(3)*(x2(1)^x(2)-x2(2)^x(2))*(x2(1)^(x(2)-1))*x(2)+0+0
> +0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0
> +0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0+0
> -dx1(2);...
> gX(2) +
> (dx1(2)*((x1(2)^(1-x(1)))/x(1))*(x1(1)^x(1)-x1(2)^x(1)
> )) +
> (dx2(2)*((x2(2)^(1-x(2)))/x(2))*(x2(1)^x(2)-x2(2)^x(2)
> )) + 0.5*x(3)*( ((x1(1)^x(1)-x1(2)^x(1))^2) +
> ((x2(1)^x(2)-x2(2)^x(2))^2) )-gX(1)];
>
> The repeated ‘+0’ is so that the each row
> of the matrix has the same number of chars.
>
> I realize looking at the code is a pain, but I would
> be very grateful for any help you can provide as I am
> really stuck.
>
> Thank you for your time,
> Ashley
Hope that's all that has to be changed.
Best wishes
Torsten.