MATLAB gets the inverse of the matrix wrong

9 views (last 30 days)
Mustafa Duran
Mustafa Duran on 24 Apr 2022
Edited: Torsten on 24 Apr 2022
Hello everyone. I write a code like this for Newton-Raphson; however system is converges. Thanks in advance to those who are interested.
clc;
clear all;
format shortG;
a=0.1;
b=0.3;
c=0.2;
d=0.2;
e=0.3;
f=0.25;
g=0.3;
h=0.2;
j=0.3;
k=0.2;
m=0.3;
n=0.2;
p=0.1;
u1=acosd(0.75);
v1=acosd(0.75);
v2=acosd(0.75);
y2=acosd(0.75);
y1=180-2*u1;
u2=180-2*u1;
max_error=0.000001;
%Initial guesses:
teta1=50;
teta2=40;
teta3=30;
teta4=60;
teta5=20;
teta6=15;
x1=0.2;
y=0.65;
for q=0:1:1
for i=1:10000
mat1=[(p+n)*cosd(teta1) e*cosd(teta2) 0 0 0 0 0 -1;
-(p+n)*sind(teta1) -e*sind(teta2) 0 0 0 0 -1 0;
0 e*cosd(teta2) 0 -c*cosd(teta1) -f*cosd(teta5) b*cosd(teta2) 0 0;
0 -e*sind(teta2) 0 c*sind(teta4) f*sind(teta5) -b*sind(teta6) 0 0;
(m-n)*cosd(teta1) 0 j*cosd(teta3) (-g-c)*cosd(teta4) -f*cosd(teta5) 0 0 0;
(-m+n)*sind(teta1) 0 -j*sind(teta3) (g+c)*sind(teta4) f*sind(teta5) 0 0 0;
k*cosd(teta1) -e*cosd(teta2) j*cosd(teta3) -h*cosd(teta4) 0 0 0 0;
-k*sind(teta1) e*sind(teta2) -j*sind(teta3) h*sind(teta4) 0 0 0 0];
f1=(p+n)*sind(teta1)+e*sind(teta2)-a*sind(q)-y;
f2=(p+n)*cosd(teta1)+e*cosd(teta2)-a*cosd(q)-x1;
f3=b*sind(teta6)-c*sind(teta4)-f*sind(teta5)+e*sind(teta2);
f4=b*cosd(teta6)-c*cosd(teta4)-f*cosd(teta5)+e*cosd(teta2);
f5=m*sind(teta1+u1)+j*sind(teta3)-g*sind(180+teta4-y2)-c*sind(teta4)-f*sind(teta5)-n*sind(teta1);
f6=m*cosd(teta1+u1)+j*cosd(teta3)-g*cosd(180+teta4-y2)-c*cosd(teta4)-f*cosd(teta5)-n*cosd(teta1);
f7=k*sind(teta1+180-y1)+j*sind(teta3)-h*sind(teta4+u2)+a*sind(q)-e*sind(teta2);
f8=k*cosd(teta1+180-y1)+j*cosd(teta3)-h*cosd(teta4+u2)+a*cosd(q)-e*cosd(teta2);
mat2=[f1;f2;f3;f4;f5;f6;f7;f8];
s=[teta1;teta2;teta3;teta4;teta5;teta6;x1;y] -inv(mat1)*mat2;
error1=abs((s(1)-teta1)/s(1));
error2=abs((s(2)-teta2)/s(2));
error3=abs((s(3)-teta3)/s(3));
error4=abs((s(4)-teta4)/s(4));
error5=abs((s(5)-teta5)/s(5));
error6=abs((s(6)-teta6)/s(6));
error7=abs((s(7)-x1)/s(7));
error8=abs((s(8)-y)/s(8));
%fprintf("mat1: %f \n",mat1(5,3))
%fprintf("s1: %f \n",s(1))
fprintf("mat1: %f \n", mat1(5,4))
fprintf("s1: %f \n",s(1))
fprintf("y: %f \n",teta1)
fprintf("error 1: %f \n",error1)
teta1=s(1);
teta2=s(2);
teta3=s(3);
teta4=s(4);
teta5=s(5);
teta6=s(6);
x1=s(7);
y=s(8);
if(error1<max_error)
if(error2<max_error)
if(error3<max_error)
if(error4<max_error)
if(error5<max_error)
if(error6<max_error)
if(error7<max_error)
if(error8<max_error)
fprintf("q=%d değeri için yardımcı koordinatlar: \n" + ...
"teta1= %f \n teta2=%f teta3= %f \n teta4=%f teta5= %f \n teta6=%f \n" + ...
"x1=%f y=%f \n ",q,teta1,teta2,teta3,teta4,teta5,teta6,x1,y)
break;
end
end
end
end
end
end
end
end
end
end
When i was searching about problem, i found that :
  • When i call to inverse of first matrix in iteration; it gives me as:
0 1.7764e-15 -4.6837 -13.863 3.7184 11.211 -6.4978 -16.025
0 0 -4.1957 -12.418 1.6467 5.415 0.61168 -1.5034
0 0 -1.7856 -5.285 0.26549 1.1085 7.6963 6.0151
0 0 -1.0181 -3.0133 1.6857 4.8475 0.23662 -1.518
0 4.4409e-16 -2.1729 -6.4312 -4.7396 -0.86534 6.4819 3.8831
1.1102e-16 -7.9502e-16 5.7564 4.1585 -5.5487 -3.5879 6.1467 4.6237
0 -1 1.8855 5.5805 -1.1721 -3.6205 1.3753 3.9725
-1 -1.1102e-16 -1.8674 -5.5271 1.0955 3.4062 -1.1124 -3.4356
When i write the matrix on another page and calculated inverse of it, results are given as:
0 1.7764e-15 -4.684 -13.863 3.7188 11.211 -6.4977 -16.025
0 8.8818e-16 -4.1956 -12.418 1.6468 5.415 0.6116 -1.5034
0 0 -1.7853 -5.2839 0.26533 1.1079 7.6962 6.0153
0 0 -1.0181 -3.0132 1.6857 4.8476 0.23658 -1.518
0 4.4409e-16 -2.1726 -6.4304 -4.7398 -0.86578 6.482 3.8834
0 8.8818e-16 5.7566 4.1589 -5.5489 -3.5882 6.1468 4.6239
0 -1 1.8855 5.5806 -1.1722 -3.6207 1.3753 3.9725
-1 4.4409e-16 -1.8675 -5.5271 1.0956 3.4064 -1.1125 -3.4357
  • As shown above, even if the value at (2,2) is 8.8818e-16, MATLAB gives 0 value. I continued about searching values, i writed 2nd matrix in iteration and take inverse of in code, it gives:
0 8.8818e-16 -4.1229 -12.51 3.1709 10.036 -5.3542 -13.878
0 1.1102e-16 -3.8266 -11.611 1.3591 5.1984 0.66694 -1.6326
0 0 -1.7972 -5.4531 0.4079 1.8427 6.8106 4.1142
0 -1.8041e-16 -0.96966 -2.9422 1.7444 4.9556 0.032418 -1.8284
0 8.8818e-16 -2.0086 -6.0945 -4.8922 -0.70095 6.241 3.0747
0 0 5.702 3.6962 -5.565 -3.3401 5.8793 3.8931
0 -1 1.7726 5.3784 -1.0455 -3.4888 1.171 3.7098
-1 3.3307e-16 -1.5747 -4.778 0.85786 2.9151 -0.78798 -2.7921
However when i write the 2nd matrix in iteration on another page and take the inverse of it, it gives:
0 0 -4.1228 -12.509 3.1707 10.035 -5.354 -13.878
0 -5.5511e-16 -3.8267 -11.611 1.3591 5.1984 0.66677 -1.6327
0 0 -1.7973 -5.4534 0.408 1.8429 6.8101 4.1136
0 1.9602e-16 -0.96963 -2.942 1.7444 4.9555 0.032273 -1.8284
0 8.8818e-16 -2.0087 -6.0949 -4.8922 -0.70068 6.2408 3.0743
0 4.4409e-16 5.702 3.6959 -5.565 -3.3399 5.8791 3.8928
0 -1 1.7726 5.3784 -1.0454 -3.4887 1.171 3.7098
-1 -4.4409e-16 -1.5747 -4.778 0.85782 2.915 -0.78798 -2.7921
  • I repeated this situation about 3rd matrix in iteration, and given inverse matrix value in iteration is:
0 -2.2204e-16 -3.8011 -11.647 2.8509 9.3197 -4.7219 -12.603
0 5.5511e-17 -3.5913 -11.004 1.148 5.0202 0.7532 -1.6835
0 -8.8818e-16 -1.7913 -5.4885 0.47991 2.277 6.3404 2.9598
0 6.9389e-18 -0.94242 -2.8876 1.8163 5.0277 -0.095593 -2.0078
0 0 -1.8857 -5.7779 -5.0559 -0.64089 6.1634 2.5847
0 0 5.7376 3.4044 -5.6734 -3.2471 5.7999 3.469
0 -1 1.7151 5.2551 -0.96483 -3.4177 1.0457 3.5597
-1 1.1102e-16 -1.3791 -4.2256 0.70029 2.5632 -0.58609 -2.3628
When i write the 3nd matrix in iteration on another page and take the inverse of it, it gives:
0 1.5543e-15 -3.801 -11.646 2.8509 9.3194 -4.7218 -12.603
-4.4409e-16 3.3404e-16 -3.5915 -11.004 1.1481 5.0203 0.75311 -1.6836
-2.2204e-16 3.058e-16 -1.7914 -5.4887 0.48003 2.2772 6.3401 2.9594
-1.1102e-16 2.3616e-16 -0.94243 -2.8875 1.8163 5.0277 -0.095579 -2.0077
-2.2204e-16 7.4989e-16 -1.8858 -5.778 -5.0558 -0.64076 6.1632 2.5844
0 0 5.7377 3.4042 -5.6735 -3.247 5.7999 3.4688
0 -1 1.7151 5.255 -0.96484 -3.4176 1.0457 3.5595
-1 -2.2204e-16 -1.3791 -4.2255 0.70031 2.5631 -0.58609 -2.3627
If i clarify openly my problem, my code is converges and i think it is about taken falsely of inverse of matrix in iteration. I think i write all of term correctly but don't know why it converges and i consider this is the problem's source.
  4 Comments
Torsten
Torsten on 24 Apr 2022
Edited: Torsten on 24 Apr 2022
Isn't
d/dx(sind(x)) = pi/180*cosd(x)
?
My advice: Never work with the degree trigonometric functions in the above context.
Mustafa Duran
Mustafa Duran on 24 Apr 2022
Thanks, I changed degree parameters to radian, however divergence problem is still exist. I considered the problem is caused at the taking inverse of matrix. Sometimes MATLAB round so small number to zero as i shown above. But i realize that divergence problem is not about it. It doesn't effect value.
I think i found the problem that is about Newton-Raphson's divergence problem:
Warning: Matrix is close to singular or badly scaled. Results
may be inaccurate. RCOND = 2.308764e-23.

Sign in to comment.

Answers (1)

Torsten
Torsten on 24 Apr 2022
Edited: Torsten on 24 Apr 2022
%Initial guesses:
teta1=50;
teta2=40;
teta3=30;
teta4=60;
teta5=20;
teta6=15;
x1=0.2;
y=0.65;
x0 = [teta1, teta2, teta3, teta4, teta5, teta6]*pi/180;
x0 = [x0,x1,y];
q = 0.1;
sol = fsolve(@(x)fun(x(1),x(2),x(3),x(4),x(5),x(6),x(7),x(8),q),x0)
function ff = fun(teta1, teta2, teta3 ,teta4 ,teta5, teta6,x1,y,q)
a=0.1;
b=0.3;
c=0.2;
d=0.2;
e=0.3;
f=0.25;
g=0.3;
h=0.2;
j=0.3;
k=0.2;
m=0.3;
n=0.2;
p=0.1;
u1=acos(0.75);
v1=acos(0.75);
v2=acos(0.75);
y2=acos(0.75);
y1=pi-2*u1;
u2=pi-2*u1;
ff(1)=(p+n)*sin(teta1)+e*sin(teta2)-a*sin(q)-y;
ff(2)=(p+n)*cos(teta1)+e*cos(teta2)-a*cos(q)-x1;
ff(3)=b*sin(teta6)-c*sin(teta4)-f*sin(teta5)+e*sin(teta2);
ff(4)=b*cos(teta6)-c*cos(teta4)-f*cos(teta5)+e*cos(teta2);
ff(5)=m*sin(teta1+u1)+j*sin(teta3)-g*sin(pi+teta4-y2)-c*sin(teta4)-f*sin(teta5)-n*sin(teta1);
ff(6)=m*cos(teta1+u1)+j*cos(teta3)-g*cos(pi+teta4-y2)-c*cos(teta4)-f*cos(teta5)-n*cos(teta1);
ff(7)=k*sin(teta1+pi-y1)+j*sin(teta3)-h*sin(teta4+u2)+a*sin(q)-e*sin(teta2);
ff(8)=k*cos(teta1+pi-y1)+j*cos(teta3)-h*cos(teta4+u2)+a*cos(q)-e*cos(teta2);
end

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!