Newtons Method but error message "Singular Matrix"
8 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Helen Hailegesesse
am 7 Dez. 2023
Bearbeitet: Torsten
am 8 Dez. 2023
When I tried to solve theta and V with Newton for systems I got this error message
- Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.765187e-24.
Previously I had no problems just solving theta with Newtons for systems, but adding V got me an ill conditioned matrix. We tried different initial values to see if they were triggering, but all of them lead to the same error message. We would be very grateful for any help we could get!
function [theta,V]= newtonfsys1(theta1,V1)
t=1; it=0; maxit=100;theta0=0;V0=200/3.6;
while norm(t)>1e-9 & it<maxit
%Creates a list of x and y values
[v1,v2,v3,v4,v5]=varden(theta0,V0);
[x0,y0]=rungekutta(v1,v2,v3,v4,v5);
[v1,v2,v3,v4,v5]=varden(theta1,V1);
[x1,y1]=rungekutta(v1,v2,v3,v4,v5);
%checking if condition is being met
for i=1:length(x1)
if (abs(11.8872-x1(i))<0.01) && (abs(18.288-x1(end))<0.01)
break
end
end
for p=1:length(x0)
if (abs(11.8872-x0(p))<0.01) && (abs(18.288-x0(end))<0.01)
break
end
end
% takes the "i" & "p" from the condition and picks a good(hopefully constant for x and y
y1=y1(i)-1.001; y0=y0(p)-1.001; x1=x1(i)-11.8872; x0=x0(p)-11.8872;
%approximate derivative
dxdtheta= (x1-x0)/(theta1-theta0);
dydtheta = (y1-y0)/(theta1-theta0);
dxdV = (x1-x0)/(V1-V0);
dydV = (y1-y0)/(V1-V0);
f=[x1;y1];
J = [dxdtheta,dxdV;dydtheta,dydV];
t=J\f;
disp(' theta V f t')
disp([theta1 norm(V1) norm(f) norm(t)])
theta0=theta1; V0=V1; theta1=theta1-t; V1=V1-t(2); it=it+1;
end
if it<maxit
theta=theta1; V=V1;
else
disp('Ingen konvergens!')
theta=[];
end
9 Kommentare
Akzeptierte Antwort
Torsten
am 7 Dez. 2023
Bearbeitet: Torsten
am 8 Dez. 2023
This is quite a tricky problem - even with the MATLAB solvers available.
For given V and theta at t = 0, integrate until x becomes 11.8872 and use the difference between y from the integration and y1 as first nonlinear equation to be satisfied. Then continue integrating until x becomes 18.288 and use the difference between y from the integration and y2 as second nonlinear equation to be satisfied.
This gives a system of two nonlinear equations in the unknowns V and theta.
I think it will become quite messy to set it up with your own ode integrator and nonlinear solver.
x0 = 0;
y0 = 2.3 ;
x1 = 11.8872;
y1 = 1.001;
x2 = 18.288;
y2 = 0;
g = 9.82;
kx = 0.01;
ky = 0.02;
m = 0.058;
fun_ode = @(t,u)[u(3);u(4);-kx*(sqrt(u(3)^2+u(4)^2))^1.5 / m;(-m*g - ky*(sqrt(u(3)^2+u(4)^2))^1.5) / m];
tspan = [0 100];
V0 = 300/3.6;
theta0 = -4*pi/180;
sol = fsolve(@(z)fun_nle(z,x0,y0,x1,y1,x2,y2,fun_ode,tspan),[V0; theta0]);
V = sol(1)
theta = sol(2)
event = @(t,u)deal([u(2),1,0]);
options = odeset('Events',event);
u0 = [x0;y0;V*cos(theta);V*sin(theta)];
[~,U] = ode45(fun_ode,tspan,u0,options);
plot(U(:,1),U(:,2))
grid on
function res = fun_nle(z,x0,y0,x1,y1,x2,y2,fun_ode,tspan)
V = z(1);
theta = z(2);
xp = V*cos(theta);
yp = V*sin(theta);
u0 = [x0;y0;xp;yp];
event1 = @(t,u)deal([u(1)-x1,1,0]);
options = odeset('Events',event1);
[~,~,TE,UE,IE] = ode45(fun_ode,tspan,u0,options);
res(1,1) = UE(2) - y1;
u0 = UE;
event2 = @(t,u)deal([u(1)-x2,1,0]);
options = odeset('Events',event2);
[~,~,TE,UE,IE] = ode45(fun_ode,tspan,u0,options);
res(2,1) = UE(2) - y2;
end
2 Kommentare
Torsten
am 7 Dez. 2023
Bearbeitet: Torsten
am 8 Dez. 2023
If you replace the calls to "fsolve" and "ode45" by calls to your own nonlinear solver and ode integrator, the code from above will also work for self-written code.
This structure - a clear separation of nonlinear solver and ode integrator - will force you not to mix parts of the ode integrator with the nonlinear solver (as you did in your code from above).
The below code e.g. replaced "fsolve" by the self-written code "nls" - simply by changing a single name:
sol = nls(@(z)fun_nle(z,x0,y0,x1,y1,x2,y2,fun_ode,tspan),[V0; theta0]);
instead of
sol = fsolve(@(z)fun_nle(z,x0,y0,x1,y1,x2,y2,fun_ode,tspan),[V0; theta0]);
Do the same for "rungekutta" instead of "ode45", and you are done.
x0 = 0;
y0 = 2.3 ;
x1 = 11.8872;
y1 = 1.001;
x2 = 18.288;
y2 = 0;
g = 9.82;
kx = 0.01;
ky = 0.02;
m = 0.058;
fun_ode = @(t,u)[u(3);u(4);-kx*(sqrt(u(3)^2+u(4)^2))^1.5 / m;(-m*g - ky*(sqrt(u(3)^2+u(4)^2))^1.5) / m];
tspan = [0 100];
V0 = 300/3.6;
theta0 = -4*pi/180;
sol = nls(@(z)fun_nle(z,x0,y0,x1,y1,x2,y2,fun_ode,tspan),[V0; theta0]);
V = sol(1)
theta = sol(2)
event = @(t,u)deal([u(2),1,0]);
options = odeset('Events',event);
u0 = [x0;y0;V*cos(theta);V*sin(theta)];
[~,U] = ode45(fun_ode,tspan,u0,options);
plot(U(:,1),U(:,2))
grid on
function res = fun_nle(z,x0,y0,x1,y1,x2,y2,fun_ode,tspan)
V = z(1);
theta = z(2);
xp = V*cos(theta);
yp = V*sin(theta);
u0 = [x0;y0;xp;yp];
event1 = @(t,u)deal([u(1)-x1,1,0]);
options = odeset('Events',event1);
[~,~,TE,UE,IE] = ode45(fun_ode,tspan,u0,options);
res(1,1) = UE(2) - y1;
u0 = UE;
event2 = @(t,u)deal([u(1)-x2,1,0]);
options = odeset('Events',event2);
[~,~,TE,UE,IE] = ode45(fun_ode,tspan,u0,options);
res(2,1) = UE(2) - y2;
end
function u = nls(f,uold)
u = zeros(numel(uold),1);
itermax = 30;
eps = 1e-6;
error = 1.0e5;
uiter = uold;
iter = 0;
while error > eps && iter < itermax
u = uiter - Jac(f,uiter)\f(uiter);
error = norm(u-uiter);
iter = iter + 1;
uiter = u;
end
end
function J = Jac(f,u)
y = f(u);
h = 1e-6;
for i = 1:numel(u)
u(i) = u(i) + h;
yh = f(u);
J(:,i) = (yh-y)/h;
u(i) = u(i) - h;
end
end
Weitere Antworten (1)
Matt J
am 7 Dez. 2023
Bearbeitet: Matt J
am 7 Dez. 2023
Landing on a point with ill-conditioned J is a hazard of plain-vanilla Newton's method. That's why people usually don't use it. They use fsolve instead.
2 Kommentare
Matt J
am 7 Dez. 2023
Bearbeitet: Matt J
am 7 Dez. 2023
Well, there are several things you can try,
- Use actual derivatives instead of finite differences
- Reduce the finite difference stepsize
- Choose a different initial point. Since it's only a 2 variable problem,you should be able to plot the objective function to see approximately where the solution lies.You could also plot cond(J) as a surface to see where the singularities are, and avoid them.
Siehe auch
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!