Error message in a program to solve a non-linear system of equations
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
Ricardo Boza Villar
am 5 Mai 2015
Kommentiert: Walter Roberson
am 8 Mai 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 Kommentare
Walter Roberson
am 6 Mai 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.
Akzeptierte Antwort
Walter Roberson
am 8 Mai 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 Kommentare
Weitere Antworten (1)
Ricardo Boza Villar
am 8 Mai 2015
2 Kommentare
Walter Roberson
am 8 Mai 2015
If that x is not close enough to [0 0 0]' then you need to lower your tol
Siehe auch
Kategorien
Mehr zu Creating and Concatenating Matrices finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!