仿Mathematica中的函数ProductLog

来源：互联网发布：crystal ball for mac 编辑：程序博客网时间：2024/04/30 13:00

今天看到Mathematica中的一个特殊函数ProductLog,即对于方程z=w*e^w中已经知道z求解w的问题，翻了下数值分析的书就用Matlab实现了。原理很简单，就是牛顿下山迭代算法，学过数值分析的应该都能写。在这里特用Matlab实现。
例子：
>> ProductLog(6+8.59*i)
ans =
1.737240895023038 + 0.618867479832110i

源代码如下：

    function wx=ProductLog(z)    if abs(real(z))+abs(imag(z))<1e-10      wx=z;    return ;    end    lnz=log(z);    zx=real(lnz);    zy=imag(lnz);    Fy=zeros(2);    F=Fy;    Fx=Fy;    temp=0;    x=zx;    y=zy;    x0=0;    y0=0;    Fy(1)=log(x^2+y^2)/2+x-zx;    Fy(2)=y-zy+atan2(y,x);    erro=abs(Fy(1))+abs(Fy(2));    loopn=1000;%最大循环次数，这里可以更改    w=1;    while loopn>0 && w>-1.05 && erro>1e-10%这里的误差限制可以更改        w=x^2+y^2;        F(1)=x/w+1;        F(2)=y/w;        w=F(1)^2+F(2)^2;        Fx(1)=(Fy(1)*F(1)-Fy(2)*F(2))/w;        Fx(2)=(Fy(2)*F(1)+Fy(1)*F(2))/w;        w=1;%下山因子        while w>=-1            x0=x-w*Fx(1);            y0=y-w*Fx(2);            Fy(1)=log(x0^2+y0^2)/2+x0-zx;            Fy(2)=y0-zy+atan2(y0,x0);            temp=abs(Fy(1))+abs(Fy(2));            if temp<erro                erro=temp;                x=x0;                y=y0;                w=-1.03;            else                 w=w-0.1;%这里的下山因子是以0.1递减的，可以更改            end        end        loopn=loopn-1;    end    wx=x+i*y;    end

0 0