牛顿拉夫森法 matlab,【原创】牛顿-拉夫森迭代求非线性方程组

function [xf,fval,iter]=newtmullt(e,x0,es,maxit,varargin)

%     written by tubehu

%     www.MATLABsky.com

%

%   input:

%     e:所求方程组

%     x0:初值,必须和求知数个数相等,否则会出错

%     es:估计误差

%     maxit:最大迭代次数

%     varargin:其它参数

%   output:

%     xf:解

%     fval:函数值

%     iter:迭代次数

%

%   example:

%   syms x1 x2

%          f=[x1^2+x1*x2-10;x1+3*x1*x2^2-57];

%          [x,f1 iter]=newtmullt(f,[1.5,3.5]',10e-6,100)

%

%   returns:

%         x =

%         1.9899

%         3.0356

%      f1 =

%        1.0e-014 *

%         -0.1776

%              0

%      iter =

%           5

%

%  or:

%           f=@(x,y,z)[x^2+x*y-z;y+3*x*y^2-57];

%           [x,f1 iter]=newtmullt(f,[1.5,3.5]',10e-6,100,10)

%

%  returns:

%        x =

%          2

%          3

%        f1 =

%          0

%          0

%        iter =

%           5

if nargin<2

error('     至少要有两个输入参数')

end

if nargin<3||isempty(es)

es=10e-6;

end

if nargin<4||isempty(maxit)

maxit=100;

end

if isa(e,'function_handle')

if nargin(e)-numel(varargin)~=numel(x0)

error('    初值必须与未知数数目相等')

end

for ii=1:nargin(e)-numel(varargin)

x{1,ii}=sym(['x',num2str(ii)]);

end

e=e(x{:},varargin{:});

clear x

end

xf=[];

x0=reshape(x0,1,[]);

xnew=x0;

iter=0;

%     try部分用于maple引擎

try

while 1

iter=iter+1;

Ja=jacobian(e,findsym(e));

J=eval(subs(Ja,findsym(e),xnew));

f=eval(subs(e,findsym(e),xnew));

xold=xnew-(J\f)';

ea=100*max(abs((J\f)./xold'));

if ea<=es||max(abs(f))<=es

xf=xold';

fval=eval(subs(e,findsym(e),xold));

break

end

xnew=xold;

if iter>maxit

break

end

end

%     catch用于mupad

catch

sy=findsym(e);

n1=findstr(sy,',');

n=numel(n1);

if isempty(n1)

x=sy;

else

x(1,1)=sym(sy(1:n1(1)-1));

for ii=2:n

x(ii,1)=sym(sy(n1(ii-1)+1:n1(ii)-1));

end

x(n+1,1)=sym(sy(n1(end)+1:end));

end

Ja=jacobian(e,x);

while 1

iter=iter+1;

J=(subs(Ja,x,xnew));

f=(subs(e,x,xnew));

xold=xnew-(J\f)';

ea=100*max(abs((J\f)./xold'));

if ea<=es||max(abs(f))<=es

fval=(subs(e,x,xold));

xf=xold';

break

end

xnew=xold;

if iter>maxit

break

end

end

end

if isempty(xf)

error('     没有找到合适的解')

end