Newtons Method but error message "Singular Matrix"

When I tried to solve theta and V with Newton for systems I got this error message

1. 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

[Torsten]

Shouldn't this be

theta1=theta1-t(1);

instead of

theta1=theta1-t;

?

[Helen Hailegesesse]

Oh thank you! That unfortunately didnt change our error.

[Torsten]

You are aware that you overwrite x0,x1,y0 and y1 as arrays by this line

% 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;

so that in the next iteration these loops

%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

are senseless because x0 and x1 have become scalars ?

Or what does rungekutta return for x0,y0,x1 and y1 ?

So you see: we are not certain how your code works. We need an executable version that reproduces the error. Otherwise it seems almost impossible to help.

[Helen Hailegesesse]

Yes we are aware, they become scalars so we can get a constant for the derivative approximation. But as you can see at the top of the loop we create a new with the help of our rungekutta (a function that does the runge kutta 4 method) x1 and y1 with our new theta1, V1 and theta0 , V0 (being the old theta1 and V1), these new x1 and y1 are not scalars therefore the condition still works. Or am I misunderstanding?I will upload the files to have an executable code where you can see the error we get. Thank you in advance!

clear all; clc; clf
format short e
theta= (-4*pi)/180; V= 300/3.6;
a(1)=(theta);a(2)=(V);


[theta,V] = newtonfsys1(theta,V)

Warning: Matrix is singular to working precision.

     theta          V              f              t
  -6.9813e-02   8.3333e+01   3.1736e+00          Inf

Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN.

     theta          V              f              t
   Inf   Inf   NaN   NaN

theta = 
   NaN

V = 
   NaN

function [theta,V]= newtonfsys1(theta1,V1)

t=1; it=0; maxit=100;theta0=1;V0=1;

while norm(t)>1e-9 & it<maxit
     %Creates a list of x and y values
    [v1,v2,v3,v4,v5]=parametrar(theta0,V0);
    [x0,y0]=rungekutta(v1,v2,v3,v4,v5);

    [v1,v2,v3,v4,v5]=parametrar(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(1); 
    V1=V1-t(2); it=it+1; 
end

if it<maxit
    theta=theta1; 
    V=V1;
else
    disp('Ingen konvergens!')
    theta=[]; 
    V=[];
end
end


function [xlista, ylista] = rungekutta(xn,yn,xpn,ypn,h) 

[xlista, ylista]= deal([]);
while yn>0
    k1 = h * ypp(VV(xpn, ypn));
    l1 = h * xpp(VV(xpn, ypn));

    k2 = h * ypp(VV(xpn + h/2, ypn + k1/2));
    l2 = h * xpp(VV(xpn + h/2, ypn + k1/2));

    k3 = h * ypp(VV(xpn + h/2, ypn + k2/2));
    l3 = h * xpp(VV(xpn + h/2, ypn + k2/2));

    k4 = h * ypp(VV(xpn + h, ypn + k3));
    l4 = h * xpp(VV(xpn + h, ypn + k3));

    k = 1/6 * (k1 + 2 .* k2 + 2 * k3 + k4);
    l = 1/6 * (l1 + 2 .* l2 + 2 * l3 + l4);

    xpn1 = xpn + l;
    ypn1 = ypn + k;

    y = yn + h  * ypn1;
    x = xn + h * xpn1;

    xlista(end+1,1)=x; ylista(end+1,1)=y;
    yn= y; xn=x; ypn=ypn1;xpn=xpn1;

end
end

function svar = ypp(V)
ky=0.02; m=0.058; g=9.82;
svar = g-(ky*V^1.5)/m;
end

function svar = xpp(V)
m=0.058; kx=0.01;
svar = -(kx*V^1.5)/m;
end

function svar = VV(xp,yp)
svar = sqrt(xp^2 + yp^2);
end

function [x0,y0,xp0,yp0,h] = parametrar(theta,V)
y0=2.30 ; 
x0=0; 
g=9.82;
kx=0.01; 
ky=0.02; 
m=0.058;
h=0.00001;
xp0=V*cos(theta); 
yp0=V*sin(theta);

end

[Helen Hailegesesse]

xpp and ypp functions are the x double dot and y double dot in the image below that describes our differential equations

[Torsten]

Sorry, but I'm completely lost in reconstructing the problem you are trying to solve with your code.

Could you try to explain ? Is it a boundary-value problem you are trying to solve with Newton's method ?

To solve the system you included, you have to rewrite it in a system of first-order differential equations:

Let

z(1) = x, z(2) = xdot, z(3) = y, z(4) = ydot.

Then your system reads

dz(1)/dt = z(2)
dz(2)/dt = -K_x*(sqrt(z(2)^2+z(4)^2))^1.5 / m
dz(3)/dt = z(4)
dz(4)/dt = (-m*g - K_y*(sqrt(z(2)^2+z(4)^2))^1.5) / m

So there should be 4 Runge-Kutta updates, not only 2 as in your code.

[Helen Hailegesesse]

Not sure if its a boundary value problem. We're trying to find the initial speed and initial angle of a tennis ball that makes it pass through the point (11.8872, 1.001) and land in (18.288, 0). Start height is y = 2.3 m. x and y describe the balls position where the x axis is ground floor, and y is vertical height. We thought to use Newtons Method for systems to solve angle and speed to fulfill those conditions.

[Helen Hailegesesse]

We arent sure of where you mean we should change our Runga-Kutta to 4 updates? Its either in our function for newton or our function for Runga-Kutta and either way we arent sure of how to implement it. Sorry for our confusion, we're very grateful for this help!

[Helen Hailegesesse]

The time of flight should be determined by the initial angle and velocity in (x0,y0). Besides ((x1,y1) = (11.8872,1.001) we also have a third point the curve needs to pass through (x2,y2) = (18.288, 0), the landing location. Its possible that there are several different combinations of angle and velocity that go through those three points, but we suspect that theres one unique curve that does. We were tasked to find the smallest angle that would achieve this.

Answer 1

[Torsten]

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]);

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

V = sol(1)

V = 5.3614e+03

theta = sol(2)

theta = -0.0793

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

[Helen Hailegesesse]

The assignment was to solve it without Matlabs integrators and solvers, that’s why we made our own Runga Kutta and Newton functions. But thank you for this extensive answer, you’ve helped a great deal!

[Torsten]

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)

V = 5.3614e+03

theta = sol(2)

theta = -0.0793

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

Answer 2

[Matt J]

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.

[Helen Hailegesesse]

Hey, thank you for answering!

We arent allowed to use fsolve, so our problem is rather to make J well conditioned using Newton's method. At the moment it stops after two iterations, and theta and V become NaN.

[Matt J]

Well, there are several things you can try,

1. Use actual derivatives instead of finite differences
2. Reduce the finite difference stepsize
3. 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.