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

1 visualizzazione (ultimi 30 giorni)
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
  5 Commenti
Walter Roberson
Walter Roberson il 6 Mag 2015
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
Ricardo Boza Villar il 6 Mag 2015
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.

Accedi per commentare.

Risposta accettata

Walter Roberson
Walter Roberson il 8 Mag 2015
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.

Più risposte (1)

Ricardo Boza Villar
Ricardo Boza Villar il 8 Mag 2015
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
  2 Commenti
Ricardo Boza Villar
Ricardo Boza Villar il 8 Mag 2015
Modificato: Ricardo Boza Villar il 8 Mag 2015
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
Walter Roberson il 8 Mag 2015
If that x is not close enough to [0 0 0]' then you need to lower your tol

Accedi per commentare.

Categorie

Scopri di più su Creating and Concatenating Matrices in Help Center e File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by