I am trying to solve 3 simultaneous nonlinear system of equations by newton's method.
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
Akshay Pratap Singh
am 23 Jan. 2019
Bearbeitet: Torsten
am 24 Jan. 2019
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 Kommentare
Akzeptierte Antwort
Torsten
am 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 Kommentare
Torsten
am 24 Jan. 2019
Bearbeitet: Torsten
am 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,:))
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Calculus finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!