Error message in a program to solve a non-linear system of equations

I would like to know why I receive an error message telling me that the matrix used to solve the problem is singular to working precision and how I could fix it.

% Function to solve a non-linear system of equations
function [x]=SENL(h)
n=3;
x=(rand(n,1)); % Random column vector
tol=10^-4;     % Relative tolerance of the solution
v=zeros(n,1);  % Zeros vector intended to create a canonical vector
v(1)=1;        %with this second line.
e=(2*tol*norm(x))*v; %This is the variable we use to solve the problem.
                     % I start it with double the value of tol*norm(v) so
                     % that the 'while' loop can start
i=0;                 % Variable 'i' for the iterations
while norm(e)>tol*norm(x) && i<50 % 2 different ways of exiting the loop
    i=i+1;
    f=fun(x); %----> It evaluates the function and returns a column vector
    J=funjac(x,h,n); %----> It calculates the jacobian using incremental
                     %      quotients. By numerical methods.
    e=J\f;           % This is used to calculate e=J^(-1)*f
    x=x-e;           % It updates the solution 'x'
    norm(e)          % norm(e) gives a smaller number with every iteration.
                     % With each iteration 'x' is closer to being the
                     % solution.
end
end

function [f]=fun(x) % It evaluates a f, which is a vector-valued function
                    % with vector-valued variables
f1=x(1)^2 +4*x(2) +3*x(3);
f2=x(1)   -x(2)^2 +x(3);
f3=9*x(1) +3*x(2) -x(3)^2;
f=[f1; f2; f3];
end

function [J]=funjac(x,h,n) % funjac returns the jacobian
                           % (differential matrix of the function)
I=ones(n);
for k=1:n
    alpha(k)=max(1,abs(x(k)))*sign(x(k));
    finc=fun(x+(alpha(k)*h*I(:,k)));
    f=fun(x);
    J(:,k)=(finc-f)/(alpha(k)*h);
end
end

[Walter Roberson]

which line is giving you that error? Is it the

e=J\f

line?

Does it give the error even on the first iteration, or does it take several iterations before it gives the error?

[Ricardo Boza Villar]

It gives the error at line 17 when I calculate:

e=J\f

The script works, but the answer is wrong, and I don't know why.

[Ricardo Boza Villar]

Can somebooty help me, please?

[Walter Roberson]

Does it give the message on the first iteration, or does it take several iterations?

As a debugging step, try displaying J and f just before doing the J\f so that they will be available on the display when the problem is encountered. And then show us a sample of the initial J and f, and show us a sample of what J and f look like when the message is generated.

[Ricardo Boza Villar]

I've made one change to the function in order to see what happens in every step. It only enters the 'while' loop once. This is the final answer:

f =

      4.6679
      0.1213
     10.0338

J =

      8.6304    8.6304    8.6304
      0.1874    0.1874    0.1874
     11.7450   11.7450   11.7450

e =

     NaN
     Inf
    -Inf

x =

     NaN
    -Inf
     Inf

i =

       1

ans =

     NaN
    -Inf
     Inf

Thanks for taking interest in my problem, I really appreciate it. See you next time.

Answer 1

[Walter Roberson]

Your initial random x is drawn from rand(), which is going to be in the range (0,1) exclusive.

Inside your funjack you have

alpha(k) = max(1,abs(x(k)) * sign(x(k));

as all of your entries are in the range 0 to 1, abs(x(k)) is going to be in the range 0 to 1, and max(1,that) is going to be 1. sign() of (0,1) is going to be 1 always, because rand() never generates exactly 0. So for the first run your alpha(k) are going to be all 1.

I cannot really trace further without knowing size(h) to know what your h*I(:,k) is intended to do. I(:,k) is going to be a column vector of 1's of length n, but your n is fixed not dependent on size(h). And since your I = ones(n) is going to be 1 everywhere, I don't understand why you bother to select the k'th column. If h is p x q then since I(:,k) is going to be n x 1, for it to work then q must be the same as n, and the result would be p x 1. Multiply that by a scalar gives p x 1, and add that to the vector "x" which is n x 1, and that is only going to work if p is 1 or p is also n. That is, the dimensions would agree if h is n x n, or if h is 1 x n. Either way you get a vector n x 1 which gets passed to fun() which will promptly use exactly 3 elements of it -- which would be bad news if n were ever changed to 2.

Skipping further analysis of the function because I am tired, if you look at your output of J you will see that all three columns are the same. And that's a singular matrix.

Answer 2

[Ricardo Boza Villar]

I don't believe it! After reading your message I realised that I had written ones(n) while in fact it was eye(n). Silly mistake.

Here's the new script and its answer to h=0.001. Thanks a lot.

% Function to solve a non-linear system of equations
function []=SENL(h)
n=3;
x=(rand(n,1)); % Random column vector
tol=10^-10;     % Relative tolerance of the solution
v=zeros(n,1);  % Zeros vector intended to create a canonical vector
v(1)=1;        %with this second line.
e=(2*tol*norm(x))*v; %This is the variable we use to solve the problem.
                     % I start it with double the value of tol*norm(v) so
                     % that the 'while' loop can start
i=0;                 % Variable 'i' for the iterations
while norm(e)>tol && i<50 % 2 different ways of exiting the loop
    i=i+1;
    f=fun(x); %----> It evaluates the function and returns a column vector
    J=funjac(x,h,n); %----> It calculates the jacobian using incremental
                     %      quotients. By numerical methods.
    e=J\f;           % This is used to calculate e=J^(-1)*f
    x=x-e;           % It updates the solution 'x'
    norm(e)          % norm(e) gives a smaller number with every iteration.
                     % With each iteration 'x' is closer to being the
                     % solution.
end
f
J
i
x
end

function [f]=fun(x) % It evaluates a f, which is a vector-valued function
                    % with vector-valued variables
f1=x(1)^2 +4*x(2) +3*x(3);
f2=x(1)   -x(2)^2 +x(3);
f3=9*x(1) +3*x(2) -x(3)^2;
f=[f1; f2; f3];
end

function [J]=funjac(x,h,n) % funjac returns the jacobian
                           % (differential matrix of the function)
I=eye(n);
for k=1:n
    alpha(k)=max(1,abs(x(k)))*sign(x(k));
    finc=fun(x+(alpha(k)*h*I(:,k)));
    f=fun(x);
    J(:,k)=(finc-f)/(alpha(k)*h);
end
    end

ans =

      1.1296

ans =

      0.2988

ans =

      0.0718

ans =

      0.0017

ans =

     1.9202e-06

ans =

     1.2262e-09

ans =

     7.8846e-13

f =

     1.0e-12 *

      0.4097
     -0.8045
     -0.8297

J =

     -0.0010    4.0000    3.0000
      1.0000   -0.0010    1.0000
      9.0000    3.0000    0.0010

i =

       7

x =

     1.0e-15 *

     -0.1675
      0.3216
     -0.3416

[Ricardo Boza Villar]

The error diminishes with every iteration.

The final f has to be [0 0 0]'

The jacobian looks fine.

The number of iterations is small but not too much.

x=[0 0 0]' is in fact a solution to the system of equations.

[Walter Roberson]

If that x is not close enough to [0 0 0]' then you need to lower your tol