Solving Nonlinear Equations Error

28 views (last 30 days)
Tim
Tim on 4 Jul 2011
I have a problem with solving two nonlinear equations using the newton method. My problem is that the code works with some values but with others I only receive the following errors:
??? Error using ==> erfc Input must be real and full.
Error in ==> normcdf at 94 p(todo) = 0.5 * erfc(-z ./ sqrt(2));
Error in ==> test4 at 14 f(1) = x*normcdf((log(x/18.5)+(0.005 +0.5*y^2))/y) - 18.5 * exp(-0.005)*normcdf((log(x/185)+(0.03 +0.5*y^2))/y)-y-0.5;
My Code is:
% Newton Raphson solution of two nonlinear algebraic equations
% set up the iteration
error1 = 1.e8;
xx(1) = 19.0; % initial guesses
xx(2) = 0.02;
iter=0;
itermax=3000.
% begin iteration
while error1>1.e-12
iter=iter+1;
x = xx(1);
y = xx(2);
% calculate the functions
f(1) = x*normcdf((log(x/18.5)+(0.005 +0.5*y^2))/y) - 18.5 * exp(-0.005)*normcdf((log(x/185)+(0.03 +0.5*y^2))/y)-y-0.5;
f(2) =(x/0.5)*normcdf((log(x/18.5)+(0.005 +0.5*y^2))/y)*y-1.42;
% calculate the Jacobian
J(1,1) = normcdf((log(x/18.5)+(0.005 +0.5*y^2))/y);
J(1,2) = x*1*(exp((-((log(x/18.5)+(0.005 +0.5*y^2))/y)^2)/2))/((sqrt(2)*pi*1)*x*y);
J(2,1) = normcdf((log(x/18.5)+(0.005 +0.5*y^2))/y)*y +(x*1*(exp((-((log(x/18.5)+(0.005 +0.5*y^2))/y)^2)/2))/((sqrt(2)*pi*1)*x*y));
J(2,2) = x*normcdf((log(x/18.5)+(0.005 +0.5*y^2))/y)+x*y*(exp((-((log(x/18.5)+(0.005 +0.5*y^2))/y)^2)/2))/((sqrt(2)*pi)*x*y)*((log(x/(18.5*exp(-0.005)))))/(y^2)+0.5;
% solve the linear equations
y = -J\f';
% move the solution, xx(k+1) - xx(k), to xx(k+1)
xx = xx + y';
% calculate norms
error1=sqrt(y(1)*y(1)+y(2)*y(2));
error(iter)=sqrt(f(1)*f(1)+f(2)*f(2));
ii(iter)=iter;
if (iter > itermax)
error1 = 0.;
s=sprintf('****Did not converge within %3.0f iterations.****',itermax);
disp(s)
end
% check if error1 < 1.e-12
end
x = xx(1);
y = xx(2);
f(1) = x*normcdf((log(x/18.5)+(0.005 +0.5*y^2))/y) - 18.5 * exp(-0.005)*normcdf((log(x/185)+(0.005 +0.5*y^2))/y)-y-0.5;
f(2) =(x/0.5)*normcdf((log(x/18.5)+(0.005 +0.5*y^2))/(y*1))*y-1.42;

Answers (0)

Categories

Find more on Symbolic Math Toolbox 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!