Newton Raphson gives answers divided by 0
2 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Jean-Paul
am 10 Mär. 2024
Bearbeitet: Matt J
am 11 Mär. 2024
I'm trying to estimate 4 angles and 2 distances using the newton raphson method with phi4 and theta as the input. I have 6 functions that equal to 0 and when i try substituting to calculate dx (the array of all my values) i get a division by 0. I transposed some arrays to columns, but id be surprised if thats the issue. The functions themselves check out and im mostly using matlab functions to make this work so I don't really know what could be wrong
Angle=3;
phi4=pi/180*Angle;
theta=linspace(0,2*pi, 72);
%First guesses
phi1=ones(1, 72);
phi2=ones(1, 72);
phi3=ones(1, 72);
phi5=ones(1, 72);
d=ones(1,72);
s=ones(1, 72);
syms phi1V phi2V phi3V phi5V dV sV;
xV=[sV, dV, phi1V, phi2V, phi3V, phi5V];
ds=0;
dd=0;
dphi1=0;
dphi2=0;
dphi3=0;
dphi4=0;
dx=[ds, dd, dphi1, dphi2, dphi3, phi5]';
lim=0.01;
i=1;
while(i<73)
%redefine functions for changes in theta
f1=25-sV+dV*cos(phi2V)+3.1*cos(phi3V);
f2=-5+dV*sin(phi2V)+3.1*sin(phi3V);
f3=-25+6.3*cos(theta(i))+27*cos(phi1V)-19.2*cos(phi2V);
f4=11+6.3*sin(theta(i))+27*sin(phi1V)-19.2*sin(phi2V);
f5=-40+6.3*cos(theta(i))+(dV-19.2)*cos(phi2V)+29.1*cos(phi5V)+16*cos(phi4);
f6=-24+6.3*sin(theta(i))+(dV-19.2)*sin(phi2V)+29.1*sin(phi5V)+16*sin(phi4);
f=[f1, f2, f3, f4, f5, f6]';
J=jacobian(f, xV);
M=-inv(J)*f;
x=[s(i), d(i), phi1(i), phi2(i), phi3(i), phi5(i)];
firstloop=1;
while(ds > lim | dd > lim | dphi1 > lim | dphi2 > lim | dphi3 > lim | dphi4 > lim | firstloop==1)
dx=subs(M, xV, x);
x=x+dx;
firstloop=0;
end
i=i+1;
end
0 Kommentare
Akzeptierte Antwort
Sam Chak
am 10 Mär. 2024
Hi @Jean-Paul
During the computation of M, certain elements involve divisions by . In the line 'dx = subs(M, xV, x)', an attempt is made to substitute the variables in 'xV = [sV, dV, phi1V, phi2V, phi3V, phi5V]' with 'x = [1, 1, 1, 1, 1, 1]', which leads to . Consequently, it triggers the 'Division-by-zero' error message.
0 Kommentare
Weitere Antworten (2)
Torsten
am 10 Mär. 2024
Bearbeitet: Torsten
am 10 Mär. 2024
Your linear systems in the course of the Newton iterations cannot be solved since the Jacobian is rank-deficient (it has only rank 5 instead of rank 6 as needed). This usually indicates that the equations you try to solve are not independent or even contradictory.
Angle=3;
phi4=pi/180*Angle;
theta=linspace(0,2*pi, 72);
%First guesses
phi1=ones(1, 72);
phi2=ones(1, 72);
phi3=ones(1, 72);
phi5=ones(1, 72);
d=ones(1,72);
s=ones(1, 72);
syms phi1V phi2V phi3V phi5V dV sV;
xV=[sV, dV, phi1V, phi2V, phi3V, phi5V];
ds=0;
dd=0;
dphi1=0;
dphi2=0;
dphi3=0;
dphi4=0;
%dx=[ds, dd, dphi1, dphi2, dphi3, phi5]';
lim=0.01;
i=1;
while(i<73)
%redefine functions for changes in theta
f1=25-sV+dV*cos(phi2V)+3.1*cos(phi3V);
f2=-5+dV*sin(phi2V)+3.1*sin(phi3V);
f3=-25+6.3*cos(theta(i))+27*cos(phi1V)-19.2*cos(phi2V);
f4=11+6.3*sin(theta(i))+27*sin(phi1V)-19.2*sin(phi2V);
f5=-40+6.3*cos(theta(i))+(dV-19.2)*cos(phi2V)+29.1*cos(phi5V)+16*cos(phi4);
f6=-24+6.3*sin(theta(i))+(dV-19.2)*sin(phi2V)+29.1*sin(phi5V)+16*sin(phi4);
f=[f1, f2, f3, f4, f5, f6]';
J=jacobian(f, xV);
%M=-inv(J)*f
x=[s(i), d(i), phi1(i), phi2(i), phi3(i), phi5(i)]
firstloop=1;
while(ds > lim | dd > lim | dphi1 > lim | dphi2 > lim | dphi3 > lim | dphi4 > lim | firstloop==1)
dx=-subs(J,xV,x)\subs(f,xV,x);%subs(M, xV, x)
rank(subs(J,xV,x))
rank([subs(J,xV,x),subs(f,xV,x)])
x=x+dx;
firstloop=0;
end
i=i+1;
end
0 Kommentare
Matt J
am 11 Mär. 2024
Bearbeitet: Matt J
am 11 Mär. 2024
Raw Newton-Raphson is a pretty terrible method, but I assume you'll tell us you have no choice...
[f,M]=preanalysis(); %Symbolic pre-analysis
%%%%Numeric stuff -- Newton-Raphson
Angle=3;
phi4=pi/180*Angle;
theta=linspace(0,2*pi, 72);
lim=1e-6;
i=1;
clear solutions
while(i<73)
F=@(x) f(x, theta(i),phi4);
Update=@(x) M(x, theta(i),phi4);
x=F(zeros(6,1))/1000; %Initial guess - needs some experimentation
dx=inf;
firstloop=1;
while( all(abs(dx)>lim) | firstloop==1)
dx=Update(x);
if any(~isfinite(dx)),
warning 'Bad update. Aborting...'
break;
end
x=x+dx;
x(1:4)=atan2( sin(x(1:4)) , cos(x(1:4))); %normalize the 2*pi periodicity
firstloop=0;
end
solutions{i}=x;
i=i+1;
end
solutions=cell2mat(solutions);
solutions(:,1:20)
function [F,M]=preanalysis()
%%%Symbolic stuff
syms phi1 phi2 phi3 phi5 d s theta phi4 real
f1=25-s+d*cos(phi2)+3.1*cos(phi3);
f2=-5+d*sin(phi2)+3.1*sin(phi3);
f3=-25+6.3*cos(theta)+27*cos(phi1)-19.2*cos(phi2);
f4=11+6.3*sin(theta)+27*sin(phi1)-19.2*sin(phi2);
f5=-40+6.3*cos(theta)+(d-19.2)*cos(phi2)+29.1*cos(phi5)+16*cos(phi4);
f6=-24+6.3*sin(theta)+(d-19.2)*sin(phi2)+29.1*sin(phi5)+16*sin(phi4);
f=[f1, f2, f3, f4, f5, f6]';
x=[phi1 phi2 phi3 phi5 d s]';
J=jacobian(f, x);
M=matlabFunction(-J\f,'Vars',{x,theta,phi4});
F=matlabFunction(f,'Vars',{x,theta,phi4});
end
0 Kommentare
Siehe auch
Kategorien
Mehr zu Stability Analysis finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!