twostep.m

偏微分方程的数值解法

function root=TwoStep(f,a,b,type,eps)
if(nargin==4)
    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));
    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)
        if(type==1)
            r1=root;
            fx1=subs(sym(f),findsym(sym(f)),r1);
            dfx=subs(sym(fun),findsym(sym(fun)),r1);
            r2=r1-fx1/dfx;
            fx2=subs(sym(f),findsym(sym(f)),r2);
            root=r2-fx2/dfx;
            tol=abs(root-r1);
        else
            r1=root;
            fx1=subs(sym(f),findsym(sym(f)),r1);
            dfx=subs(sym(fun),findsym(sym(fun)),r1);
            r2=r1-fx1/dfx;
            fx2=subs(sym(f),findsym(sym(f)),r2);
            root=r2-fx2*(r2-r1)/(2*fx2-fx1);
            tol=abs(root-r1);
        end
    end
end