Error message in a program to solve a non-linear system of equations
1 view (last 30 days)
Show older comments
Ricardo Boza Villar
on 5 May 2015
Commented: Walter Roberson
on 8 May 2015
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 Comments
Walter Roberson
on 6 May 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.
Accepted Answer
Walter Roberson
on 8 May 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.
0 Comments
More Answers (1)
Ricardo Boza Villar
on 8 May 2015
2 Comments
Walter Roberson
on 8 May 2015
If that x is not close enough to [0 0 0]' then you need to lower your tol
See Also
Categories
Find more on Dynamic System Models in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!