今天看到 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-10wx=z;return ;endlnz=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>=-1x0=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<erroerro=temp;x=x0;y=y0;w=-1.03;elsew=w-0.1;%這里的下山因子是以0.1遞減的，可以更改endendloopn=loopn-1;endwx=x+i*y;end

總結

以上是生活随笔為你收集整理的仿Mathematica中的函数ProductLog的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。