I am trying to solve 3 simultaneous nonlinear system of equations by newton's method.
3 views (last 30 days)
Show older comments
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
0 Comments
Accepted Answer
Torsten
on 23 Jan 2019
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 Comments
Torsten
on 24 Jan 2019
Edited: Torsten
on 24 Jan 2019
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,:))
More Answers (0)
See Also
Categories
Find more on Calculus in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!