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

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

Answer 1

[Torsten]

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

[Akshay Pratap Singh]

"jacobian" is coming ok. there is no problem with it.

[Torsten]

And what happens to "a" ? "J" and "g" have to be evaluated at "a" in each iteration of your while-loop. That's why I recommend to do all calculations in your while-loop numerically, not symbolically.

[Torsten]

This might work (I didn't test it):

... {Code as provided}
g = [e1; e2; e3];
J = jacobian([e1, e2, e3], [x, y, z]);
while del > 1e-5

   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

[Akshay Pratap Singh]

Thanks Torsten, it worked.

[Akshay Pratap Singh]

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]

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,:))