基于matlab的非线性薛定谔方程的数值算法研究

­1­基于MATLAB的非线性薛定谔方程的数值算法研究雷超四川农业大学农业工程系，四川雅安（625014）E­mail（[e-mail]）摘要：本文简单介绍了一种基于MATLAB来求解非线性薛定谔方程的分步傅立叶数值算法。非线性薛定谔方程是由线性算符和非线性算符两部分组成，依据分步傅立叶数值算法的思想，当光场每通过一段微小距离的时候，先发生线性算符代表的色散数学方程，再计算非线性算符代表的非线性效应数学方程，从而得到近似的数值结果。在计算存在高阶偏微分项时，在时域中计算非常不便，可利用傅立叶变换，把偏微分方程变换为代数方程来运算，进而求出其相应的数值结果。从算法的编写角度出发，该算法在MATLAB中易于实现，而且比较简单清晰，能够广泛的应用在非线性薛定谔方程的数值求解之中。关键词：非线性薛定谔方程；傅立叶变换；分步傅立叶数值算法；MATLAB中图分类号：O4151.引言非线性薛定谔方程是研究光脉冲在光纤中传输的基本方程，大多数文献都是直接引用非线性薛定谔方程，采用分步傅里叶法数值求解非线性薛定谔方程。本文从通过MATLAB编程来求解非线性薛定谔方程，并给出了分步傅里叶法数值求解算法程序。MATLAB是美国MathWorks公司自20世纪80年代中期推出的数学软件，优秀的数值计算能力和卓越的数据可视化能力使其很快在数学软件中脱颖而出。随着版本的不断升级，它在数值计算及符号计算功能上得到了进一步完善。到目前为止，MATLAB已经发展成为多学科、多种工作平台的功能强大的大型软件。2.非线性薛定谔方程非线性薛定谔方程，简称NLS方程，是一个非线性偏微分方程，通常情况下无法求出解析解，只能求出它的数值解。在含有非线性色散介质地脉冲传输问题中应用非常广泛。非线性薛定谔方程的基本形式为:（2­1）uuuiuxxt2||2其中是未知的复值函数.u3.分步傅立叶数值算法目前，采用分步傅立叶算法(SplitstepFourierMethod)求解非线性薛定谔方程的数值解应用比较多。分步傅立叶方法最早是在1937年开始应用的，这种方法己经被证明是相同精度下数值求解非线性薛定愕方程最快的方法，部分原因是它采用了快速傅立叶变换算法(FastFourierTransformAlgorithm)。基于MATLAB科学计算软件以及MATLAB强大的符号计算功能，完全可以实现分步傅立叶数值算法来对脉冲形状和频谱进行仿真。­2­一般情况下，光脉冲沿光纤传播时受到色散和非线性效应的共同作用，假设当传输距离很小的时候，两者相互独立作用，那么，根据这种思想可建立如下分步傅立叶数值算法的数学模型：把待求解的非线性薛定谔方程写成以下形式:（3­1）UNDzU)ˆˆ(其中是线性算符，代表介质的色散和损耗，是非线性算符，它决定了脉冲传输过程中DˆNˆ光纤的非线性效应。一般来讲，沿光纤的长度方向，色散和非线性是同时作用的。分步傅立叶法假设在传输过程中，光场每通过一小段距离，色散和非线性效应可以分别作用，得到近似结果。也就h是说脉冲从到的传输过程中分两步进行。第一步，只有非线性作用，方程(2)式中的zhz；第二步，再考虑线性作用，方程(2)式中的.0ˆD0ˆN这样方程(2)在这两步中可分别简化为:（3­2）UNzUUDzUˆˆ得到了上面两个方程（3­2），就可以分别求解非线性作用方程和线性作用方程，然后讨论分步傅立叶法的数值算法。由于方程（3­2）是一个偏微分方程，需要通过傅立叶变换把偏微分方程转换为代数方程，进行运算。傅立叶变换的定义如下：（3­3）dTTizUTzUzUFdTTiTzUzUTzUF)exp(),(~21),()],(~[)exp(),(),(~)],([1在计算时一般采用快速傅立叶变换（FFT）算。为了保证精度要求，一般还需要)],([TzUF反复调整纵向传输步长z和横向脉冲取样点数T来保证计算精度。4.分步傅立叶数值算法的MATLAB实现为了研究分步傅立叶法在Matlab上的实现，下面通过举例来说明Matlab在求解的非线性薛定谔方程的优越性。现待求解的非线性薛定谔方程如下:（4­1）0||42222AAiTAiAzA其中，是光场慢变复振幅，z是脉冲沿光纤传播的距离；),(TzA,是群速度；是色散系数；是非线性系数；gvztT/1,11gv)/(kmps)/1(kmw)是光纤损耗系数，它与用分贝表示的损耗系数的关系为：)/1(km)/(kmdBdB.343.4dB首先，可以将方程（4­1）归一化振幅：，是入射脉冲的峰值功率，0/),(PTzAU0P此时方程(4­1)可改写为:­3­（4­2）UUiPTUiUzU2022||42为了使用分步傅立叶法求解方程(4­2)，将方程(4­2)写成以下形式:UNDzU)ˆˆ(进一步，可以得出如下方程（4­3）：（4­3）2022||ˆ22ˆUiPNTUiD然后，按照步骤step1和步骤step2，依次计算方程（4­3）的线性算符和非线性算符。最后在步骤step3中，运行步骤step1和步骤step2的MATLAB程序，得出线性算符和非线性算符的精确数值解及其仿真曲线。Step1线性算符方程的求解线性算符的方程如下：Dˆ（4­4）UTUizU2222用傅立叶变换方法，得到一个常微分方程（4­5）：（4­5）UiiUzU~4)(~2~2解方程(4­5)得:（4­6）]42exp[),0(~),(~2ziUzU式中是初值的傅立叶变换，将进行反傅立叶变换就得到了),0(~U),0(TU),(~zU。方程(4­6)的求解公式为:),(TzU（4­7）)]},0([]2)2{exp[(~),(2TUFziFTzU其中和分别表示傅立叶变换和反傅立叶变换运算。FF~Step2非线性算符方程的求解非线性部分的方程如下：Nˆ（4­8）UUiPzU20||同Step1的方法，解方程（4­8），得到：（4­9）]|),0(|exp[),0(~),(~20zTUiPUzU式中是初值的傅立叶变换，将进行反傅立叶变换就得到了),0(~U),0(TU),(~zU。方程(4­9)的求解公式为:),(TzU（4­10）)]},0([]|),0(|{exp[~),(20TUFzTUiPFTzU­4­其中和分别表示傅立叶变换和反傅立叶变换运算。FF~Step3算法在MATLAB中的实现在Matlab中，设有限时长序列的长度为，它对应于一个频域内的长)(nx)1(NnN度为N的有限长序列.对应的角频，其中)1)((NkkX)(kXNTkk2)()1(NkT是序列的采样时间间隔.)(nx这种正反离散傅立叶变换的关系式为:（4­11）)1()2exp()(1)]([)()1()2exp()()]([)(11NnnkNjkXNkXIDFTnxNknkNjnxnxDFTkXNjNj然后用Matlab中的离散傅立叶变换(DFT)函数fft和离散傅立叶反变换(IDFT)的函数ifft来实现方程(4­4)，(4­8)式中的傅立叶和反傅立叶变换运算。进一步，得到方程(4­7)，(4­10)的数值解及仿真曲线。实现分步傅立叶变换算法MATLAB的脚本文件程序如下：Po；%输入光强，单位Walpha；%光纤损耗值，单位dB/kmgamma；%光纤非线性参数to；%初始脉冲宽度，单位秒C；%第一次计算输入的啁啾参数b2；%波数的倒数cputime=0;tic;ln=1;i=sqrt(­1);pi=3.1415926535;alph=alpha/(4.343);Ld=(to^2)/(abs(b2));%扩散长度，单位是mAo=sqrt(Po);%光振幅tau=­4096e­12:1e­12:4095e­12;dt=1e­12;h=1000;%步长forii=0.1:0.1:1.5%不同的光纤长度不同，这个量可变z=ii*Ld;u=Ao*exp(­((1+i*(­C))/2)*(tau/to).^2);figure(1)plot(abs(u),'r');title('InputPulse');xlabel('Time');ylabel('Amplitude');gridon;holdon;l=max(size(u));­5­具体在实现分步傅立叶变换算法的时候，除了上述脚本本件定义的参数外，还可以根据实际情况，修改程序部分参数，直到达到你想要求解精度和曲线。最后，通过测试一组参数，得到方程（4­1）在该算法下的MATLAB运算结果。MATLAB总共用时177.2660s，求得的的结果曲线如下图1所示。结果表明，算法正确而且精度也比较高，能够在非线性薛定谔方程的求解中广泛应用。fwhm1=find(abs(u)abs(max(u)/2));fwhm1=length(fwhm1);dw=1/l/dt*2*pi;w=(­1*l/2:1:l/2­1)*dw;u=fftshift(u);%零延迟对中的谱w=fftshift(w);%零延迟对中的谱spectrum=fft(fftshift(u));%快速离散傅立叶变换forjj=h:h:zspectrum=spectrum.*exp(g1);%g1为线性算符e的指数表达式f=ifft(spectrum);%快速离散反傅立叶变换f=f.*exp(g2);%g2为非线性算符e的指数表达式spectrum=fft(f);%快速离散傅立叶变换spectrum=spectrum.*exp(g1);endf=ifft(spectrum);%快速离散反傅立叶变换op_pulse(ln,:)=abs(f);%保存在所有间隔点上的输出脉冲fwhm=find(abs(f)abs(max(f)/2));fwhm=length(fwhm);ratio=fwhm/fwhm1;pbratio(ln)=ratio;dd=atand((abs(imag(f)))/(abs(real(f))));phadisp(ln)=dd;%保存脉冲相位ln=ln+1;endtoc;cputime=toc;figure(2);mesh(op_pulse(1:1:ln­1,:));title('PulseEvolution');xlabel('Time');ylabel('distance');zlabel('amplitude');figure(3)plot(pbratio(1:1:ln­1),'k');xlabel('Numberofsteps');ylabel('Pulsebroadeningratio');gridon;holdon;figure(5)plot(phadisp(1:1:ln­1),'k');xlabel('distancetravelled');ylabel('phasechange');gridon;holdon;disp('CPUtime:'),disp(cputime);­6­图1.仿真结果参考文献[1]程雪苹.直接微扰方法和非线性薛定谔方程族的近似解[D].浙江：浙江师范大学,2006.[2]齐志勇.掺铒光纤放大器的模拟器研究[D].湖北：华中科技大学,2004.[3]施娟.基于对称分步傅立叶算法的光孤子仿真［J］.技术前沿,2008,第10卷（第1期）[4]张兴坊,郑义,李爱萍,徐云峰.皮秒脉冲在光子晶体光纤中的压缩效应［J］.激光技术,200,第31卷（第3