非线性薛定谔方程数值解的MATLAB仿真

[键入作者姓名][键入文档标题]——利用分步快速傅里叶变换对光纤中光信号的传输方程进行数值求解1、非线性薛定谔方程非线性薛定谔方程(nonlinear Schrodinger equation ，NLSE)是奥地利物理学家薛定谔于1926 年提出的，应用在量子力学系统中。

由于量子力学主要研究粒子的动力学运动状态，所以不能运用牛顿力学公式来表示。

通常在量子力学中，研究系统的状态一般通过波函数(x ，t)来表示。

而对波函数的研究主要是求解非线性薛定谔方程。

本文主要研究光脉冲在光纤中传输状态下的演变。

一般情况下，光脉冲信号在光纤中传输时，同时受到光纤的色散和非线性效应的影响。

通过Maxwell 方程，考虑到光纤的色散和非线性效应，可以推导出光信号在光纤中的传输方程，即非线性薛定谔方程。

NLSE 是非线性偏微分方程，一般很难直接求出解析解，于是通过数值方法进行求解。

具体分为两大类：(1)分布有限差分法(split-step finite differencemethod ，SSFD)；(2)分步傅里叶变换法(split-step Fourier transform method ，SSFT)。

一般情况，在达到相同精度，由于分步傅里叶变换法采用运算速度快的快速傅里叶变换，所以相比较有限差分法运算速度快一到两个数量级。

于是本文介绍分步傅里叶变换法来对光纤中光信号的传输方程，即非线性薛定谔方程进行数值求解。

并通过MATLAB 软件对结果数值仿真。

非线性薛定谔方程的基本形式为:22||t xx iu u u u =+其中u 是未知的复值函数.目前，采用分步傅立叶算法(Split step Fourier Method)求解非线性薛定谔方程的数值解应用比较多。

分步傅立叶方法最早是在1937年开始应用的，这种方法己经被证明是相同精度下数值求解非线性薛定愕方程最快的方法，部分原因是它采用了快速傅立叶变换算法(Fast Fourier Transform Algorithm)。

基于MATLAB 科学计算软件以及MATLAB 强大的符号计算功能，完全可以实现分步傅立叶数值算法来对脉冲形状和频谱进行仿真。

一般情况下，光脉冲沿光纤传播时受到色散和非线性效应的共同作用，假设当传输距离很小的时候，两者相互独立作用，那么，根据这种思想可建立如下分步傅立叶数值算法的数 学模型：（I ）把待求解的非线性薛定谔方程写成以下形式:ˆˆ()U D N U z∂=+∂ 其中ˆD是线性算符，代表介质的色散和损耗， ˆN 是非线性算符，它决定了脉冲传输过程中光纤的非线性效应。

一般来讲，沿光纤的长度方向，色散和非线性是同时作用的。

分步傅立叶法假设在传输过程中，光场每通过一小段距离h ，色散和非线性效应可以分别作用，得到近似结果。

也就是说脉冲从z 到z +h 的传输过程中分两步进行。

第一步，只有非线性作用，方程(II)式中的ˆD=0；第二步，再考虑线性作用，方程(II)式中的ˆN =0 这样方程(2)在这两步中可分别简化为:ˆˆU DU zU N U z∂=⋅∂∂=⋅∂得到了上面两个方程（III ），就可以分别求解非线性作用方程和线性作用方程，然后讨论分步傅立叶法的数值算法。

由于方程（III ）是一个偏微分方程，需要通过傅立叶变换把偏微分方程转换为代数方 程，进行运算。

傅立叶变换的定义如下：1[(,)](,)(,)exp()1[(,)](,)(,)exp()2F U z T U z U z T i T dTF U z U z T U z i T dT ωωωωωπ+∞-∞+∞--∞⎧⎪==⎪⎨⎪==-⎪⎩⎰⎰ 在计算[(,)]F U z T 时一般采用快速傅立叶变换（FFT ）算。

为了保证精度要求，一般还需要反复调整纵向传输步长z 和横向脉冲取样点数T 来保证计算精度。

2、分步傅立叶数值算法的MATLAB 实现现待求解的非线性薛定谔方程如下:222024A i A A i A A z T αβγ∂∂+--=∂∂ 其中，A (z ,T )是光场慢变复振幅，z 是脉冲沿光纤传播的距离；（II ）（III ）（IV ） （V ）1T t z β=-,11/g v β=,v g 是群速度；(/)ps km β是色散系数；(1/)w km γ⋅是非线性系数；(1/)km α是光纤损耗系数，它与用分贝表示的损耗系数(/)dB dB km α的关系为：4.343dB αα=.首先，可以将方程（V ）归一化振幅：0(,)/U A z T P =， 0P 是入射脉冲的峰值功率，此时方程(V)可改写为：220224U i U U P i U U z Tαβγ∂∂=-++∂∂ 为了使用分步傅立叶法求解方程(VI )，将方程(VI )写成以下形式:ˆˆ()U DN U z∂=+∂ 进一步，可以得出如下方程（VII ）：2222ˆ2ˆi UT DN P i U βαγ∂-+∂==然后，按照步骤1和步骤2，依次计算方程（VII ）的线性算符和非线性算符。

最后在步骤3 中，运行步骤1和步骤2的MATLAB 程序，得出线性算符和非线性算符的精确数值解及其仿真曲线。

步骤1 线性算符方程的求解线性算符的方程如下：2222i UUT U zβα∂-+∂∂=∂用傅立叶变换方法，得到一个常微分方程（IX ）：2()24U i i U U z αωβ∂=--∂ 解方程(IX)得:22(,)(0,)exp[]4i U z U z βωαωω-=式中(0,)U ω是初值(0,)U T 的傅立叶变换，将(,)U z ω进行反傅立叶变换就得到了（VI ）（VII ）（VIII ）（IX ）（X ）(,)U z T 。

方程(X)的求解公式为:2(,){exp[()][(0,)]}22i zU z T F F U T βωα=-⋅其中F 和F 分别表示傅立叶变换和反傅立叶变换运算。

步骤2 非线性算符方程的求解非线性部分的方程如下：20U P i U U zγ∂=∂ 同Step1的方法，解方程（XII ），得到：20(,)(0,)exp[(0,)]U z U Pi U T z ωωγ=式中(0,)U ω是初值(0,)U T 的傅立叶变换，将(,)U z ω进行反傅立叶变换就得到了(,)U z T 。

方程(XIII)的求解公式为:20(,){exp[(0,)][(0,)]}U z T F Pi U T z F U T γ=⋅其中F 和F 分别表示傅立叶变换和反傅立叶变换运算。

步骤3 算法在MATLAB 中的实现在Matlab 中，设有限时长序列()x n 的长度为(1)N n N ≤≤，它对应于一个频域内的长度为N 的有限长序列()(1)X k k N ≤≤，()X k 的角频2()(1)kk k N NTπω=≤≤，其中T 是序列()x n 的采样时间间隔.这种正反离散傅立叶变换的关系式为:112()[()]()exp()(1)12()[()]()exp()(1)Nj Nj X k DFT x n x n j k n k N Nx n IDFT X k X k j k n n N NNππ====-⋅⋅⋅ ≤≤==⋅⋅⋅ ≤≤∑∑ 然后用Matlab 中的离散傅立叶变换(DFT )函数fft 和离散傅立叶反变换(IDFT )的函数ifft 来实现方程(VIII)，(XII)式中的傅立叶和反傅立叶变换运算。

进一步，得到方程(XI)，(XIV) 的数值解及仿真曲线。

最后，通过测试一组参数，得到方程（V ）在该算法下的MATLAB 运算结果。

（XI ）（XII ）（XIII ）（XV ）（XIV ）MATLAB总共用时34.26s，求得的的结果曲线如下图所示。

结果表明，算法正确而且精度也比较高，能够在非线性薛定谔方程的求解中广泛应用。

附录MATLAB的脚本文件源代码：Po=200; %输入光强，单位Walpha=3.5; %光纤损耗值，单位dB/kmgamma=20; %光纤非线性参数to=1; %初始脉冲宽度，单位秒C=50; %第一次计算输入的啁啾参数b2=1000; %波数的倒数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;%步长for ii=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('Input Pulse'); xlabel('Time'); ylabel('Amplitude');grid on;hold on;l=max(size(u));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)); %快速离散傅立叶变换for jj=h:h:zspectrum=spectrum.*exp(g1); %g1为线性算符e的指数表达式f=ifft(spectrum); %快速离散反傅立叶变换f=f.*exp(g2);%g2为非线性算符e的指数表达式spectrum=fft(f); %快速离散傅立叶变换f=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('Pulse Evolution');xlabel('Time'); ylabel('distance'); zlabel('amplitude'); figure(3)plot(pbratio(1:1:ln-1),'k');xlabel('Number of steps');ylabel('Pulse broadening ratio');grid on;hold on;figure(4)plot(phadisp(1:1:ln-1),'k');xlabel('distance travelled');ylabel('phase change');grid on;hold on;disp('CPU time:'), disp(cputime);。