Windows LiveWindows Live™
 
 
impreza's profileLostPhotosBlogListsMore Tools Help

Blog
    April 20

    非线性方程求根

                         
     
    二分法:          
    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)
        n=n+1;
        a1=a2;
        r1=root;
        root=subs(sym(f),findsym(sym(f)),r1)+r1;
        x(n)=root;
        if(n>2)
            m=m+1;
            a2=x(m)-(x(m+1)-x(m))^2/(x(m+2)-2*x(m+1)+x(m));
            tol=abs(a2-a1);
        end  
    end
    root=a2;
    >> [r,n]=AtkenStablePoint('1/sqrt(x)+x-2',0.5)
     
    r =
     
        0.3820
     
     
    n =
     
         4
    Steffensen加速迭代：
    function [root,n]=StevenStablePoint(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;
        y=subs(sym(f),findsym(sym(f)),r1)+r1;
        z=subs(sym(f),findsym(sym(f)),y)+y;
        root=r1-(y-r1)^2/(z-2*y+r1);
        tol=abs(root-r1);
    end
    >> [r,n]=StevenStablePoint('1/sqrt(x)+x-2',0.5)
     
    r =
     
        0.3820
     
     
    n =
     
         4
    弦截法：       
    function root=Secant(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;
        fa=subs(sym(f),findsym(sym(f)),a);
        fb=subs(sym(f),findsym(sym(f)),b); 
        root=a-(b-a)*fa/(fb-fa);
        while(tol>eps)
            r1=root;
            fx=subs(sym(f),findsym(sym(f)),r1);
            s=fx*fa;
            if(s==0)
                root=r1;
            else
                if(s>0)
                    root=b-(r1-b)*fb/(fx-fb);
                else
                    root=a-(r1-a)*fa/(fx-fa);
                end
            end
            tol=abs(root-r1);
        end
    end
    >> Secant('x^3-3*x+1',0,1)
     
    ans =
     
        0.3473
    Steffensen弦截法：
    function root=StevenSecant(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;
        fa=subs(sym(f),findsym(sym(f)),a);
        fb=subs(sym(f),findsym(sym(f)),b);
        faa=subs(sym(f),findsym(sym(f)),fa+a);
        root=a-fa*fa/(faa-fa);
        while(tol>eps)
            r1=root;
            fx=subs(sym(f),findsym(sym(f)),r1);
            v=fx+r1;
            fxx=subs(sym(f),findsym(sym(f)),v);
            root=r1-fx*fx/(fxx-fx);
            tol=abs(root-r1);
        end
    end
    >> StevenSecant('x^3-3*x+1',0,1)
     
    ans =
     
        0.3473
    抛物线法：
    function root=Parabola(f,a,b,x,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;
        fa=subs(sym(f),findsym(sym(f)),a);
        fb=subs(sym(f),findsym(sym(f)),b);
        fx=subs(sym(f),findsym(sym(f)),x);
        d1=(fb-fa)/(b-a);
        d2=(fx-fb)/(x-b);
        d3=(d2-d1)/(x-a);
        B=d2+d3*(x-b);
        root=x-2*fx/(B+sign(B)*sqrt(B^2-4*fx*d3));
        t=zeros(3);
        t(1)=a;
        t(2)=b;
        t(3)=x;
        while(tol>eps)
            t(1)=t(2);
            t(2)=t(3);
            t(3)=root;
            f1=subs(sym(f),findsym(sym(f)),t(1));
            f2=subs(sym(f),findsym(sym(f)),t(2));
            f3=subs(sym(f),findsym(sym(f)),t(3));
            d1=(f2-f1)/(t(2)-t(1));
            d2=(f3-f2)/(t(3)-t(2));
            d3=(d2-d1)/(t(3)-t(1));
            B=d2+d3*(t(3)-t(2));
            root=t(3)-2*f3/(B+sign(B)*sqrt(B^2-4*f3*d3)); 
            tol=abs(root-t(3));
        end
    end
    >> Parabola('sqrt(x)+log(x)-2',1,4,2)
     
    ans =
     
        1.8773
    Newton法：
    function root=NewtonRoot(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));
        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-fa/dfa;
        else
            root=b-fb/dfb;
        end
        while(tol>eps)
            r1=root;
            fx=subs(sym(f),findsym(sym(f)),r1);
            dfx=subs(sym(fun),findsym(sym(fun)),r1);
            root=r1-fx/dfx;
            tol=abs(root-r1);
        end
    end
    >> NewtonRoot('sqrt(x)-x^3+2',0.5,2)
     
    ans =
     
        1.4759
    简化Newton法：
    function root=SimpleNewton(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));
        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)
            df0=1/dfa;
            root=a-df0*fa;
        else
            df0=1/dfb;
            root=b-df0*fb;
        end
        while(tol>eps)
            r1=root;  
            fx=subs(sym(f),findsym(sym(f)),r1);
            root=r1-df0*fx;
            tol=abs(root-r1);
        end
    end
    >> SimpleNewton('sqrt(x)-x^3+2',1.2,2)
     
    ans =
     
        1.4759
    Newton下山法：
    function root=NewtonDown(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));
        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);
            toldf=1;
            alpha=2;
            while toldf>0
                alpha=alpha/2;
                root=r1-alpha*fx/dfx;
                fv=subs(sym(f),findsym(sym(f)),root);
                toldf=abs(fv)-abs(fx);
            end
            tol=abs(root-r1);
        end
    end
    >> NewtonDown('sqrt(x)-x^3+2',1.2,2)
     
    ans =
     
    1.4759
     
    两步迭代法：
    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
    >> Two_Step('log(x)-sin(x)+1',0.1,3,1)
     
    ans =
     
        0.7013
    重根的迭代法：
    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
    >> MultiRoot('(sin(x)-x+2)*x*x',-2,3,1e-8)
     
    ans =
     
         0
     
    非线性方程组的数值解法
     
    不动点迭代：
    function [r,n]=mulStablePoint(myf,x0,eps)
    if nargin==1
        eps=1.0e-4;
    end
     
    r=myf(x0);
    n=1;
    tol=1;
    while tol>eps
        x0=r;
        r=myf(x0);
        tol=norm(r-x0);
        n=n+1;
        if(n>100000)
            disp('迭代步数太多，可能不收敛！');
            return;
        end
    end
    Newton法：
    function [r,n]=mulNewton(myf,dmyf,x0,eps)
    if nargin==1
        eps=1.0e-4;
    end
     
    r=x0-myf(x0)/dmyf(x0);
    n=1;
    tol=1;
    while tol>eps
        x0=r;
        r=x0-myf(x0)/dmyf(x0);
        tol=norm(r-x0);
        n=n+1;
        if(n>100000)
            disp('迭代步数太多，可能不收敛！');
            return;
        end
    end
    简化Newton法：
    function [r,n]=mulSimNewton(myf,dmyf,x0,eps)
    if nargin==1
        eps=1.0e-4;
    end
     
    r=x0-myf(x0)/dmyf(x0);
    c=dmyf(x0);
    n=1;
    tol=1;
    while tol>eps
        x0=r;
        r=x0-myf(x0)/c;
        tol=norm(r-x0);
        n=n+1;
        if(n>100000)
            disp('迭代步数太多，可能不收敛！');
            return;
        end
    end
    Newton下山法：
    function [r,n]=mulDNewton(myf,dmyf,x0,eps)
    if nargin==1
        eps=1.0e-4;
    end
     
    r=x0-myf(x0)/dmyf(x0);
    n=1;
    tol=1;
    while tol>eps
        x0=r;
        ttol=1;
        w=1;
        F1=norm(myf(x0));
        while ttol>=0
            r=x0-w*myf(x0)/dmyf(x0);
            ttol=norm(myf(r))-F1;
            w=w/2;
        end
        tol=norm(r-x0);
        n=n+1;
        if(n>100000)
            disp('迭代步数太多，可能不收敛！');
            return;
        end
    end
    拟Newton法：
    function [r,n]=mulVNewton(myf,x0,A,eps)
    if nargin==1 
        A=eye(length(x0));
    else
        if nargin==2
            eps=1.0e-4;
        end
    end
     
    r=x0-myf(x0)/A;
    n=1;
    tol=1;
    while tol>eps
        x0=r;
        r=x0-myf(x0)/A;
        y=r-x0;
        z=myf(r)-myf(x0);
        A1=A+(z-y*A)'*y/norm(y);
        A=A1;
        n=n+1;
        if(n>100000)
            disp('迭代步数太多，可能不收敛！');
            return;
        end
        tol=norm(r-x0);
    end                                       
    12:55 PM | Add a comment | Permalink | Blog it
    April 13

    总算有点眉目了

         当一个动态的协同节点上有某存储的消息与查询消息匹配时，生成的响应消息要在传统意义上原路返回几乎不可能，所以可以采取基于梯度场的机制。梯度值的如何计算是一个问题（跳数、延迟应该是主要考虑因素）。此外，有多个查询消息匹配同一个目标消息时，如何优化传输？若一个查询得到多个响应，如何过滤？整个过程的模型如何建立呢？
         先找全所有变量再说吧！HOHOHOOH。。

     

     

    12:00 PM | Add a comment | Permalink | Blog it