syms r1 r2 r3 r4 m1 m2 Z;f1=r2; f1r1=diff(f1,r1);f1r2=diff(f1,r2);f1r3=diff(f1,r3);f1r4=diff(f1,r4);f1m1=diff(f1,m1);f1m2=diff(f1,m2) %Calculation of partial derivatives
f2=r3;f2r1=diff(f2,r1);f2r2=diff(f2,r2);f2r3=diff(f2,r3);f2r4=diff(f2,r4);f2m1=diff(f2,m1);f2m2=diff(f2,m2) %Calculation of partial derivativesf3=r4f3r1=diff(f3,r1);f3r2=diff(f3,r2);f3r3=diff(f3,r3);f3r4=diff(f3,r4);f3m1=diff(f3,m1);f3m2=diff(f3,m2) %Calculation of partial derivativesf4=Z^2*r3-(Z^4-2*Z^3+Z^2)*m2f4r1=diff(f4,r1);f4r2=diff(f4,r2);f4r3=diff(f4,r3);f4r4=diff(f4,r4);f4m1=diff(f4,m1);f4m2=diff(f4,m2) %Calculation of partial derivativesf5=m2;f5r1=diff(f5,r1);f5r2=diff(f5,r2);f5r3=diff(f5,r3);f5r4=diff(f5,r4);f5m1=diff(f5,m1);f5m2=diff(f5,m2) %Calculation of partial derivativesf6=m1*r2+(Z/2-1)*r1*m2;f6r1=diff(f6,r1);f6r2=diff(f6,r2);f6r3=diff(f6,r3);f6r4=diff(f6,r4);f6m1=diff(f6,m1);f6m2=diff(f6,m2) %Calculation of partial derivatives% n=input('Enter the number of decimal places:');
epsilon = 1*10^-(5)err=2r10=0;r20=0;r30=rand-10; %rand
r40=rand-10; %randm10=1;m20=rand-10; %randr2inf=Z^2;r3inf=0;m1inf=0;r1k(1)=r10; r2k(1)=r20; r3k(1)=r30; r4k(1)=r40; m1k(1)=m10; m2k(1)=m20; %The intial approximations
while (err)>epsilon r10=0;r20=0;r30=r30+0.0001; %randr40=r40+0.0001; %randm10=1;m20=m20+0.0001; %randr1k(1)=r10; r2k(1)=r20; r3k(1)=r30; r4k(1)=r40; m1k(1)=m10; m2k(1)=m20; %The intial approximationsfor i=1:100 f1k=vpa(subs(f1,[r1 r2 r3 r4 m1 m2 ],[r1k(i) r2k(i) r3k(i) r4k(i) m1k(i) m2k(i)])); % Calculation of function values at each iteration
f2k=vpa(subs(f2,[r1 r2 r3 r4 m1 m2 ],[r1k(i) r2k(i) r3k(i) r4k(i) m1k(i) m2k(i)])); % Calculation of function values at each iterationf3k=vpa(subs(f3,[r1 r2 r3 r4 m1 m2 ],[r1k(i) r2k(i) r3k(i) r4k(i) m1k(i) m2k(i)])); % Calculation of function values at each iterationf4k=vpa(subs(f4,[r1 r2 r3 r4 m1 m2 ],[r1k(i) r2k(i) r3k(i) r4k(i) m1k(i) m2k(i)])); % Calculation of function values at each iterationf5k=vpa(subs(f5,[r1 r2 r3 r4 m1 m2 ],[r1k(i) r2k(i) r3k(i) r4k(i) m1k(i) m2k(i)])); % Calculation of function values at each iterationf6k=vpa(subs(f6,[r1 r2 r3 r4 m1 m2 ],[r1k(i) r2k(i) r3k(i) r4k(i) m1k(i) m2k(i)])); % Calculation of function values at each iterationf1r1k=vpa(subs(f1r1,[r1 r2 r3 r4 m1 m2 ],[r1k(i) r2k(i) r3k(i) r4k(i) m1k(i) m2k(i)])); f1r2k=vpa(subs(f1r2,[r1 r2 r3 r4 m1 m2 ],[r1k(i) r2k(i) r3k(i) r4k(i) m1k(i) m2k(i)])); f1r3k=vpa(subs(f1r3,[r1 r2 r3 r4 m1 m2 ],[r1k(i) r2k(i) r3k(i) r4k(i) m1k(i) m2k(i)]));f1r4k=vpa(subs(f1r4,[r1 r2 r3 r4 m1 m2 ],[r1k(i) r2k(i) r3k(i) r4k(i) m1k(i) m2k(i)]));f1m1k=vpa(subs(f1m1,[r1 r2 r3 r4 m1 m2 ],[r1k(i) r2k(i) r3k(i) r4k(i) m1k(i) m2k(i)]));f1m2k=vpa(subs(f1m2,[r1 r2 r3 r4 m1 m2 ],[r1k(i) r2k(i) r3k(i) r4k(i) m1k(i) m2k(i)]));f2r1k=vpa(subs(f2r1,[r1 r2 r3 r4 m1 m2 ],[r1k(i) r2k(i) r3k(i) r4k(i) m1k(i) m2k(i)])); f2r2k=vpa(subs(f2r2,[r1 r2 r3 r4 m1 m2 ],[r1k(i) r2k(i) r3k(i) r4k(i) m1k(i) m2k(i)])); f2r3k=vpa(subs(f2r3,[r1 r2 r3 r4 m1 m2 ],[r1k(i) r2k(i) r3k(i) r4k(i) m1k(i) m2k(i)]));f2r4k=vpa(subs(f2r4,[r1 r2 r3 r4 m1 m2 ],[r1k(i) r2k(i) r3k(i) r4k(i) m1k(i) m2k(i)]));f2m1k=vpa(subs(f2m1,[r1 r2 r3 r4 m1 m2 ],[r1k(i) r2k(i) r3k(i) r4k(i) m1k(i) m2k(i)]));f2m2k=vpa(subs(f2m2,[r1 r2 r3 r4 m1 m2 ],[r1k(i) r2k(i) r3k(i) r4k(i) m1k(i) m2k(i)]));f3r1k=vpa(subs(f3r1,[r1 r2 r3 r4 m1 m2 ],[r1k(i) r2k(i) r3k(i) r4k(i) m1k(i) m2k(i)])); ;f3r2k=vpa(subs(f3r2,[r1 r2 r3 r4 m1 m2 ],[r1k(i) r2k(i) r3k(i) r4k(i) m1k(i) m2k(i)])); f3r3k=vpa(subs(f3r3,[r1 r2 r3 r4 m1 m2 ],[r1k(i) r2k(i) r3k(i) r4k(i) m1k(i) m2k(i)]));f3r4k=vpa(subs(f3r4,[r1 r2 r3 r4 m1 m2 ],[r1k(i) r2k(i) r3k(i) r4k(i) m1k(i) m2k(i)]));f3m1k=vpa(subs(f3m1,[r1 r2 r3 r4 m1 m2 ],[r1k(i) r2k(i) r3k(i) r4k(i) m1k(i) m2k(i)]));f3m2k=vpa(subs(f3m2,[r1 r2 r3 r4 m1 m2 ],[r1k(i) r2k(i) r3k(i) r4k(i) m1k(i) m2k(i)]));f4r1k=vpa(subs(f4r1,[r1 r2 r3 r4 m1 m2 ],[r1k(i) r2k(i) r3k(i) r4k(i) m1k(i) m2k(i)])); f4r2k=vpa(subs(f4r2,[r1 r2 r3 r4 m1 m2 ],[r1k(i) r2k(i) r3k(i) r4k(i) m1k(i) m2k(i)])); f4r3k=vpa(subs(f4r3,[r1 r2 r3 r4 m1 m2 ],[r1k(i) r2k(i) r3k(i) r4k(i) m1k(i) m2k(i)]));f4r4k=vpa(subs(f4r4,[r1 r2 r3 r4 m1 m2 ],[r1k(i) r2k(i) r3k(i) r4k(i) m1k(i) m2k(i)]));f4m1k=vpa(subs(f4m1,[r1 r2 r3 r4 m1 m2 ],[r1k(i) r2k(i) r3k(i) r4k(i) m1k(i) m2k(i)]));f4m2k=vpa(subs(f4m2,[r1 r2 r3 r4 m1 m2 ],[r1k(i) r2k(i) r3k(i) r4k(i) m1k(i) m2k(i)]));f5r1k=vpa(subs(f5r1,[r1 r2 r3 r4 m1 m2 ],[r1k(i) r2k(i) r3k(i) r4k(i) m1k(i) m2k(i)])); f5r2k=vpa(subs(f5r2,[r1 r2 r3 r4 m1 m2 ],[r1k(i) r2k(i) r3k(i) r4k(i) m1k(i) m2k(i)])); f5r3k=vpa(subs(f5r3,[r1 r2 r3 r4 m1 m2 ],[r1k(i) r2k(i) r3k(i) r4k(i) m1k(i) m2k(i)]));f5r4k=vpa(subs(f5r4,[r1 r2 r3 r4 m1 m2 ],[r1k(i) r2k(i) r3k(i) r4k(i) m1k(i) m2k(i)]));f5m1k=vpa(subs(f5m1,[r1 r2 r3 r4 m1 m2 ],[r1k(i) r2k(i) r3k(i) r4k(i) m1k(i) m2k(i)]));f5m2k=vpa(subs(f5m2,[r1 r2 r3 r4 m1 m2 ],[r1k(i) r2k(i) r3k(i) r4k(i) m1k(i) m2k(i)]));f6r1k=vpa(subs(f6r1,[r1 r2 r3 r4 m1 m2 ],[r1k(i) r2k(i) r3k(i) r4k(i) m1k(i) m2k(i)])); f6r2k=vpa(subs(f6r2,[r1 r2 r3 r4 m1 m2 ],[r1k(i) r2k(i) r3k(i) r4k(i) m1k(i) m2k(i)])); f6r3k=vpa(subs(f6r3,[r1 r2 r3 r4 m1 m2 ],[r1k(i) r2k(i) r3k(i) r4k(i) m1k(i) m2k(i)]));f6r4k=vpa(subs(f6r4,[r1 r2 r3 r4 m1 m2 ],[r1k(i) r2k(i) r3k(i) r4k(i) m1k(i) m2k(i)]));f6m1k=vpa(subs(f6m1,[r1 r2 r3 r4 m1 m2 ],[r1k(i) r2k(i) r3k(i) r4k(i) m1k(i) m2k(i)]));f6m2k=vpa(subs(f6m2,[r1 r2 r3 r4 m1 m2 ],[r1k(i) r2k(i) r3k(i) r4k(i) m1k(i) m2k(i)]));xyk = [r1k r2k r3k r4k m1k m2k ]'; %The old values of x,y,z
J = [f1r1k f1r2k f1r3k f1r4k f1m1k f1m2k; f2r1k f2r2k f2r3k f2r4k f2m1k f2m2k; f3r1k f3r2k f3r3k f3r4k f3m1k f3m2k; f4r1k f4r2k f4r3k f4r4k f4m1k f4m2k; f5r1k f5r2k f5r3k f5r4k f5m1k f5m2k; f6r1k f6r2k f6r3k f6r4k f6m1k f6m2k]; % The jacobian matrix
Fk = [f1k f2k f3k f4k f5k f6k]'; % Matrix of function values
MATLAB: Why inverse of Jacobian matrix FAIL?
jacobinnewton raphsonnonlinearnumericalode
Best Answer