I am trying to solve 3 simultaneous nonlinear system of equations by newton's method.

1 次查看（过去 30 天）
Akshay Pratap Singh
编辑： Torsten 2019-1-24
I am trying to solve 3 simultaneous nonlinear system of equations by newton's method. I am getting an error in solving that system.
The code is written below.
clear all
syms x y z
a=[1;1;1];
% The Newton-Raphson iterations starts here
del=1;
indx=0;
h=4;
gma=18.4;
ka1=0.2;
kp1=8.76;
ka2=0.2;
kp2=8.76;
sind=0.437;
cosd=0.9;
pa1=ka1*gma*(h*x+0.5*(x^2));
pp1=kp1*gma*0.5*(x^2);
pa2=ka2*gma*(x*y+0.5*(y^2));
pp2=kp2*gma*(y*(h+x)+0.5*(y^2));
za1=(0.5*h*(x^2)+(x^3)/6)/(h*x+0.5*(x^2));
zp2=(0.5*(h+x)*(y^2)+((y^3)/3))/((h+x)*y+0.5*(y^2));
za2=(0.5*x*(y^2)+((y^3)/3))+(x*y+0.5*(y^2));
e1=(pp1*sind)-(pa1*sind)-(pp2*sind)-(pa2*sind);
e2=pp1*cosd+pa2*cosd-pa1*cosd-pp2*cosd-z;
e3=pp1*cosd*(x/3)+pp2*cosd*zp2-pa1*cosd*za1-pa2*cosd*za2-z*(x+(h/3));
while del>1e-6
g=[e1; e2; e3];
J=jacobian([e1, e2, e3], [x, y, z]);
delx=-inv(J)*g;
a=a+delx;
del=max(abs(g));
indx=indx+1;
end
'NEWTON-RAPHSON SOLUTION CONVERGES IN ITERATIONS',indx,pause
'FINAL VALUES OF x ARE',x

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2019-1-23
Don't use symbolic variables together with Newton-Raphson.
Instead of using "jacobian" , look at the solution suggested by John D'Errico under
https://de.mathworks.com/matlabcentral/answers/28066-numerical-jacobian-in-matlab
  6 个评论
显示 4更早的评论
Akshay Pratap Singh
Akshay Pratap Singh 2019-1-24
Thank you very much Torsten
Now, I am having one problem that the code is working fine upto h=12 and giving good results. but as I increase h to 13, 14, 15 and so on. It is giving wrong values. Your help will be highly appreciated. Thank you in advance.
The code is written below as well as attached-
clear all
syms x y z
a=[1;1;1];
format longEng
% The Newton-Raphson iterations starts here
del=1;
indx=0;
h=13;
gma=18.4;
ka1=0.2;
kp1=8.7;
ka2=0.32;
kp2=8.7;
sinda1=0.437;
sindp1=0.437;
sinda2=0.437;
sindp2=0.437;
cosda1=0.9;
cosdp1=0.9;
cosda2=0.9;
cosdp2=0.9;
pa1=ka1*gma*0.5*(x^2);
pp1=kp1*gma*0.5*(x^2);
pa2=ka2*gma*((x*y)+0.5*(y^2));
pp2=kp2*gma*((x*y)+0.5*(y^2));
za1=x/3;
zp1=x/3;
zp2=(0.5*x*(y^2)+((y^3)/3))/((x*y)+0.5*(y^2));
za2=(0.5*x*(y^2)+((y^3)/3))/((x*y)+0.5*(y^2));
e1=(pp1*sindp1)-(pa1*sinda1)-(pp2*sindp2)-(pa2*sinda2);
e2=(pp1*cosdp1)+(pa2*cosda2)-(pa1*cosda1)-(pp2*cosdp2)-z;
e3=(pp1*cosdp1*zp1)+(pp2*cosdp2*zp2)-(pa1*cosda1*za1)-(pa2*cosda2*za2)-(z*(x+h));
g=[e1; e2; e3];
J=jacobian([e1, e2, e3], [x, y, z]);
while del>1e-6
gnum = double(subs(g,[x,y,z],[a(1),a(2),a(3)]));
Jnum = double(subs(J,[x,y,z],[a(1),a(2),a(3)]));
delx = -Jnum\gnum;
a = a + delx;
del = max(abs(gnum));
indx = indx + 1;
end
'NEWTON-RAPHSON SOLUTION CONVERGES IN ITERATIONS',indx,pause
'FINAL VALUES OF a ARE',a
Torsten
Torsten 2019-1-24
编辑：Torsten 2019-1-24
Use the result from the previous h as initial guess for the next h.
In your case:
Use
a=[1.87;0.74;17.47]; % approximate result for h=12
as initial "a" for h=13.
And I think you get faster convergence if you replace
gnum = double(subs(g,[x,y,z],[a(1),a(2),a(3)]));
Jnum = double(subs(J,[x,y,z],[a(1),a(2),a(3)]));
by
gnum = vpa(subs(g,[x,y,z],[a(1),a(2),a(3)]));
Jnum = vpa(subs(J,[x,y,z],[a(1),a(2),a(3)]));
You could try this code:
syms x y z h
a=[1;1;1];
format longEng
% The Newton-Raphson iterations starts here
H=linspace(2,30,29);
gma=18.4;
ka1=0.2;
kp1=8.7;
ka2=0.32;
kp2=8.7;
sinda1=0.437;
sindp1=0.437;
sinda2=0.437;
sindp2=0.437;
cosda1=0.9;
cosdp1=0.9;
cosda2=0.9;
cosdp2=0.9;
pa1=ka1*gma*0.5*(x^2);
pp1=kp1*gma*0.5*(x^2);
pa2=ka2*gma*((x*y)+0.5*(y^2));
pp2=kp2*gma*((x*y)+0.5*(y^2));
za1=x/3;
zp1=x/3;
zp2=(0.5*x*(y^2)+((y^3)/3))/((x*y)+0.5*(y^2));
za2=(0.5*x*(y^2)+((y^3)/3))/((x*y)+0.5*(y^2));
e1=(pp1*sindp1)-(pa1*sinda1)-(pp2*sindp2)-(pa2*sinda2);
e2=(pp1*cosdp1)+(pa2*cosda2)-(pa1*cosda1)-(pp2*cosdp2)-z;
e3=(pp1*cosdp1*zp1)+(pp2*cosdp2*zp2)-(pa1*cosda1*za1)-(pa2*cosda2*za2)-(z*(x+h));
g=[e1; e2; e3];
J=jacobian([e1, e2, e3], [x, y, z]);
A=zeros(3,numel(H));
for i=1:numel(H)
del = 1.0;
indx = 0;
while del>1e-6
gnum = vpa(subs(g,[x,y,z,h],[a(1),a(2),a(3),H(i)]));
Jnum = vpa(subs(J,[x,y,z,h],[a(1),a(2),a(3),H(i)]));
delx = -Jnum\gnum;
a = a + delx;
del = max(abs(gnum));
indx = indx + 1;
end
A(:,i)=double(a)
end
plot(H,A(3,:))

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Systems of Nonlinear Equations 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by