why inverse of Jacobian matrix FAIL??

18 views (last 30 days)
jannat alsaidi on 16 Nov 2019
Commented: Walter Roberson on 18 Nov 2019
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 derivatives
f3=r4
f3r1=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 derivatives
f4=Z^2*r3-(Z^4-2*Z^3+Z^2)*m2
f4r1=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 derivatives
f5=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 derivatives
f6=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=2
r10=0;
r20=0;
r30=rand-10; %rand
r40=rand-10; %rand
m10=1;
m20=rand-10; %rand
r2inf=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; %rand
r40=r40+0.0001; %rand
m10=1;
m20=m20+0.0001; %rand
r1k(1)=r10; r2k(1)=r20; r3k(1)=r30; r4k(1)=r40; m1k(1)=m10; m2k(1)=m20; %The intial approximations
for 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 iteration
f3k=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 iteration
f4k=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 iteration
f5k=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 iteration
f6k=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 iteration
f1r1k=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
jannat alsaidi on 16 Nov 2019
I shared a part of the code, here is the code that needed for jacobian,
syms r1 r2 r3 r4 m1 m2 ;
Z=0;
for w=1:11
f=[r2; r3;r4;Z^2*r3-(Z^4-2*Z^3+Z^2)*m2;m2;m1*r2+(Z/2-1)*r1*m2];
for i=1:6
J(i,1)=diff(f(i),r1);
J(i,2)=diff(f(i),r2);
J(i,3)=diff(f(i),r3);
J(i,4)=diff(f(i),r4);
J(i,5)=diff(f(i),m1);
J(i,6)=diff(f(i),m2);
end
Z=Z+0.1;
end
J
inv(J)

Walter Roberson on 16 Nov 2019
J =
[ 0, 1, 0, 0, 0, 0]
[ 0, 0, 1, 0, 0, 0]
[ 0, 0, 0, 1, 0, 0]
[ 0, 0, 1, 0, 0, 0]
[ 0, 0, 0, 0, 0, 1]
[ -m2/2, m1, 0, 0, r2, -r1/2]
Notice that the second and fourth lines are the same. J does not have full rank.
Walter Roberson on 18 Nov 2019
Please check my transcription. In particular check the variable the functions are in, and check the naming of the squiggle . Does the squiggle happen to be the variable the function is over, f(xi) instead of f(t) ?
syms f(t) theta(t) xi
df = diff(f); d2f = diff(df); d3f = diff(d2f); d4f = diff(d3f);
dtheta = diff(theta); d2theta = diff(dtheta);
eqn1 = d4f == xi^2 * d2f - (xi^4 - 2*xi^3 + xi^2)*dtheta
eqn2 = d2theta == theta * df + (xi/2 - 1) * f * dtheta;