Matlab非线性方程求根

 

                                                        非线性方程求根                     

 

二分法 :           

function root=MultiRoot(f,a,b,eps)

if(nargin==3)

    eps=1.0e-4;

end

 

f1=subs(sym(f),findsym(sym(f)),a);

f2=subs(sym(f),findsym(sym(f)),b);

if(f1==0)

    root=a;

end

if(f2==0)

    root=b;

end

 

if(f1*f2>0)

    disp(' 两端点函数值乘积大于0!');

    return;

else

    tol=1;

    fun=diff(sym(f));

    ddf=diff(sym(fun));

    fa=subs(sym(f),findsym(sym(f)),a);

    fb=subs(sym(f),findsym(sym(f)),b);

    dfa=subs(sym(fun),findsym(sym(fun)),a);

    dfb=subs(sym(fun),findsym(sym(fun)),b);

    if(dfa>dfb)

        root=a;

    else

        root=b;

    end

    while(tol>eps)

        r1=root;

        fx=subs(sym(f),findsym(sym(f)),r1);

        dfx=subs(sym(fun),findsym(sym(fun)),r1);

        ddfx=subs(sym(ddf),findsym(sym(ddf)),r1);

        root=r1-fx*dfx/(dfx*dfx-fx*ddfx);

        tol=abs(root-r1);

    end

end

Example:

>> HalfInterval('x^3-3*x+1',0,1)

 

ans =

 

    0.3473

黄金分割法 :

function root=hj(f,a,b,eps)

if(nargin==3)

    eps=1.0e-4;

end

 

f1=subs(sym(f),findsym(sym(f)),a);

f2=subs(sym(f),findsym(sym(f)),b);

if(f1==0)

    root=a;

end

if(f2==0)

    root=b;

end

 

if(f1*f2>0)

    disp(' 两端点函数值乘积大于0!');

    return;

else   

    t1=a+(b-a)*0.382;

    t2=a+(b-a)*0.618;

    f_1=subs(sym(f),findsym(sym(f)),t1);

    f_2=subs(sym(f),findsym(sym(f)),t2);

    tol=abs(t1-t2);

    while(tol>eps)                           % 精度控制

        if(f_1*f_2<0)

            a=t1;

            b=t2;

        else

            fa=subs(sym(f),findsym(sym(f)),a);        

            if(f_1*fa>0)

                a=t2;

            else

                b=t1;

            end

        end

        t1=a+(b-a)*0.382;

        t2=a+(b-a)*0.618;

        f_1=subs(sym(f),findsym(sym(f)),t1);

        f_2=subs(sym(f),findsym(sym(f)),t2);

        tol=abs(t2-t1);

    end

    root=(t1+t2)/2;                             % 输出根

end

 

>> hj('x^3-3*x+1',0,1)

 

ans =

 

    0.3473

不动点迭代:

function [root,n]=StablePoint(f,x0,eps)

if(nargin==2)

    eps=1.0e-4;

end

 

tol=1;

root=x0;

n=0;

while(tol>eps)

    n=n+1;

    r1=root;

    root=subs(sym(f),findsym(sym(f)),r1)+r1;

    tol=abs(root-r1);

end

>> [r,n]=StablePoint('1/sqrt(x)+x-2',0.5)

 

r =

 

    0.3820

 

 

n =

 

     4

Atken 加速迭代法：

function [root,n]=AtkenStablePoint(f,x0,eps)

if(nargin==2)

    eps=1.0e-4;

end

 

tol=1;

root=x0;

x(1:2)=0;

n=0;

m=0;

a2=x0;

while(tol>eps)