why inverse of Jacobian matrix FAIL??

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

[Walter Roberson]

What you posted is missing two end statements.

What you posted uses r1k(i) in a for loop with i potentially being > 1, and it defines r1k(1) but it never assigns to any other r1k() location, so r1k(2) does not exist.

[jannat alsaidi]

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)

Answer 1

[Walter Roberson]

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.

[jannat alsaidi]

What shall I do to have a jacobian matrix with rank =6, these functions (f) came from a reduction of fourth and second order nonlinear system of ODEs, if you have an idea, help plz.

[Walter Roberson]

What are the original equations?

[jannat alsaidi]

[Walter Roberson]

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;