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

1、admin非线性薛定谓方程数值解的 MATLAB仿真利用分步快速傅里叶变换对光纤中光信号的传输方程进行数值求解1、非线性薛定调方程非线性薛定谓方程(nonlinear Schrodinger equation , NLSE谑奥地利物理学家薛定渭 于1926年提出的，应用在量子力学系统中。由于量子力学主要研究粒子的动力学运动状态，所以不能运用牛顿力学公式来表示。通常在量子力学中，研究系统的状态一般通过波函数(x,t)来表示。而对波函数的研究主要是求解非线性薛定谓方程。本文主要研究光脉冲在光纤中传输状态下的演变。一般情况下，光脉冲信号在光纤中传输时，同时受到光纤的色散和非线性效应的影响。 通过Ma

2、xwell方程，考虑到光纤的色散和非线性效应，可以推导出光信号在光纤中的传输 方程，即非线性薛定谓方程。NLSE是非线性偏微分方程，一般很难直接求出解析解，于是通过数值方法进行求解。具体分为两大类：(1)分布有限差分法(split-step finitedifferencemethod , SSFD); (2)分步傅里口t变换法 (split-step Fourier transform method, SSFT)。一般情况，在达到相同精度，由于分步傅里叶变换法采用运算速度快的快速傅里叶 变换，所以相比较有限差分法运算速度快一到两个数量级。于是本文介绍分步傅里叶变换法来对光纤中光信号的传输方程

3、，即非线性薛定谓方程进行数值求解。并通过MATLA啾件对结果数值仿真。非线性薛定谓方程的基本形式为：iut =uxx +2| u |2 u(I)其中u是未知的复值函数.目前，采用分步傅立叶算法 (Split step Fourier Method)求解非线性薛定谓方程的数值解应用比较多。分步傅立叶方法最早是在 1937年开始应用的，这种方法己经被证明是相同 精度下数值求解非线性薛定愕方程最快的方法，部分原因是它采用了快速傅立叶变换算法 (FastFourier Transform Algorithm) 。基于 MATLAB学计算软件以及 MATLAB1 大的符号计 算功能，完全可以实现分步傅立

4、叶数值算法来对脉冲形状和频谱进行仿真。一般情况下，光脉冲沿光纤传播时受到色散和非线性效应的共同作用，假设当传输距离很小的时候，两者相互独立作用，那么，根据这种思想可建立如下分步傅立叶数值算法的数 学模型：把待求解的非线性薛定谓方程写成以下形式：电=(3 +的U(II)二 zN?是非线性算符，它决定了脉冲传输过其中D?是线性算符，代表介质的色散和损耗, 程中光纤的非线性效应。一般来讲，沿光纤的长度方向，色散和非线性是同时作用的。分步傅立叶法假设在传输过程中，光场每通过一小段距离h,色散和非线性效应可以分别作用，得到近似结果。也就是说脉冲从z到z+h的传输过程中分两步进行。第一步，只有非线性作用，

5、方程 （II）式中的8二0；第二步，再考虑线性作用，方程 （ii）式中的N? =o这样方程（2）在这两步中可分别简化为：U=D? U(III)；z=N? U；z得到了上面两个方程 （III）,就可以分别求解非线性作用方程和线性作用方程，然后讨论分步傅立叶法的数值算法。由于方程（III）是一个偏微分方程，需要通过傅立叶变换把偏微分方程转换为代数方 程，进行运算。傅立叶变换的定义如下:*-boFU(z,T) =U(z, ) = U(z,T)exp(i T)dT(IV)FU (z, ) =U(z,T) =,_'.U (z, )exp(-i T)dT在计算FU（z,T）时一般采用快速傅立叶变换

6、（FFT）算。为了保证精度要求，一般还 需要反复调整纵向传输步长 z和横向脉冲取样点数T 来保证计算精度。2、分步傅立叶数值算法的MATLA限现现待求解的非线性薛定谓方程如下 ：4 A.:z 2-i；：2A4；:T222'i A A = 0(V)其中，A（z,T）是光场慢变复振幅，z是脉冲沿光纤传播的距离；T = t -P1z,P1 =1/Vg ,vg是群速度；P（ps/km）是色散系数；丫（1/ w km）是非线性系数；口（1/ km）是光纤损耗系数，它与用分贝表示的损耗系数udB（dB/km）的关系为：c（dB = 4.343t.首先，可以将方程（V）归一化振幅：U = A（z,T

7、）/Jp", P0是入射脉冲的峰值功率,此时方程（V）可改写为:U:z一 U2U24 FT2+ ¥P0iUU(VI)为了使用分步傅立叶法求解方程（VI）,将方程（VI）写成以下形式U =（D? + N?）U；z进一步，可以得出如下方程（ VII ）: i ；：2U- -I ''2D?二2-1-2N? = VF0i U 2(VII )然后，按照步骤1和步骤2,依次计算方程（VII ）的线性算符和非线性算符。最后在步 骤3中，运行步骤1和步骤2的MATLA皿序，得出线性算符和非线性算符的精确数值解及其 仿真曲线。步骤1线性算符方程的求解线性算符的方程如下:.:U

8、.:za +-i ；：2U2 2-Hu 2(VIII )用傅立叶变换方法，得到一个常微分方程（ IX）:U:z二 i(i -)2 ：U - U24(IX)解方程（IX）得:(X)_ i ；，2 -2：U (z, ) =U (0, )expz4式中U （0,0）是初值U （0,T）的傅立叶变换，将 U （z,c）进行反傅立叶变换就得到了U（z,T）。方程（X）的求解公式为：U（z,T） = Fexp（-归2 -«）- FU （0,T）（XI）22其中F和F分别表示傅立叶变换和反傅立叶变换运算。步骤2非线性算符方程的求解非线性部分的方程如下:rU2，、= YP0i U U(XII)；z同

9、Stepl的方法，解方程(XII),得到：E2U(z2) =U(0,co)exp¥P0i U(0,T) z(XIII)式中U(0, 是初值U (0,T)的傅立叶变换，将U(z,co)进行反傅立叶变换就得到了 U (z,T)。方程(XIII)的求解公式为： E 一一2一U (z,T) =Fexp ¥P0i U (0,T) z FU (0,T)(XIV)其中F和F分别表示傅立叶变换和反傅立叶变换运算。步骤3算法在MATLAB的实现在Matlab中，设有限时长序列 x(n)的长度为N(1 <n < N),它对应于一个频域内的长度2二 k为N的有限长序列X(k)(1Wk

10、MN), X(k)的角频®(k)=(1 Ek E N),其中T是序列NTx(n)的采样时间间隔.这种正反离散傅立叶变换的关系式为NX(k) = DFTx(n) = " x(n)exp(-jj=1k n) N(1 < k < N)1 J2 二x(n) =IDFTX(k) X(k)exp( j 一 k n) (1 < n < N)N jdN(XV)然后用Matlab中的离散傅立叶变换(DFT)函数fft和离散傅立叶反变换(IDFT)的函数ifft来实现方程(VIII) , (XII)式中的傅立叶和反傅立叶变换运算。进一步，得到方程(XI) , (XIV)

11、的数值解及仿真曲线。最后，通过测试一组参数， 得到方程(V)在该算法下的MATLAB运算结果。MATLAB总共用时34.26 s,求得的的结果曲线如下图所示。结果表明，算法正确而且精度也比较高，能够在非线性薛定谓方程的求解中广泛应用。6 4o.o.8E.4,21.8 1111o.2一兵BUcBPGOq Bs-nd0.2 L 0510Number of steps15Pulse Evolution0.050.04、0.03 .0.02.-0.0101610400010000180006000distance2000Time50distance travelled45 40 .36 -30 -26

12、 -20 15 -附录MATLA的脚本文件源代码: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=(toA2)/(abs(b2); % 扩散长度，单位是 mAo=sqrt(Po); % 光振幅tau=-4096e-12:1e-12:4095e-12;dt=1e-12;h=1000;% 步

13、长for ii=0.1:0.1:1.5 %不同的光纤长度不同，这个量可变 z=ii*Ld;u=Ao*exp(-(1+i*(-C)/2)*(tau/to)A2);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=(-

14、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); %快速离散傅立叶变换 spectrum=spectrum.*exp(g1);endf=ifft(spectrum); %快速离散反傅立

15、叶变换 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