muller.m

来自「余向阳博士Matlab教程的源代码。 包括：1.非线性方程的解法；2.线性方程」· M 代码 · 共 51 行

function [p,Y,err]=muller(f1043,p0,p1,p2,delta,max1)
%  f是给定的非线性方程。
%  p0,p1,p2是给定的初始值。
%  delta是近似值的误差界。
%  max1是迭代次数的上限。
%  p为所求得到方程的近似解。
%  y=f(p)。
%  err为p的误差估计。
P=[p0,p1,p2],
k=0,
Y=[feval('f1043',p0),feval('f1043',p1),feval('f1043',p2)],
err=abs(p2-p1);
for k=1:max1
   h0=P(1)-P(3);
   h1=P(2)-P(1);
   e0=Y(1)-Y(3);
   e1=Y(2)-Y(3);
   c=Y(3);
   denom=h1*h0^2-h0*h1^2;
   a=(e0*h1-e1*h0)/denom;
   b=(e1*h0^2-e0*h1^2)/denom;
   if b^2-4*a*c>0
        disc=sqrt(b^2-4*a*c);
   else 
           disc=0;
end
if b<0
          disc=-disc;
    end
  z=-2*c/(b+disc);
  p=P(3)+z;
  if abs(p-P(2))< abs(p-P(1)) 
Q=[P(2) P(1) P(3)];
      P=Q;
            Y=[feval('f1043',P(1)),feval('f1043',P(2)),feval('f1043',P(3))];
  end
  if abs(p-P(3))<abs(p-P(2))
            R=[P(1)  P(3)  P(2)];
      P=R;
            Y=[feval('f1043',P(1)) ,feval('f1043',P(2)) ,feval('f1043',P(3)) ];
  end
  P(3)=p;
  Y(3)=feval('f1043',P(3));
  y=Y(3);
  err=abs(z);
    p,err,k,y,
  if(err<delta)|(abs(Y)==0)
           break
        end
end