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

基于 MATLAB 的非线性薛定谔方程的数值算法研究 如果您需要更多论文可以到www.docin.com/week114进行免费查阅。 基于 MATLAB 的非线性薛定谔方程的数值算法研究 雷超 四川农业大学农业系，四川雅安 (625014) E-mail(leichaocool123@sina.com) 摘 要:本文简单介绍了一种基于 MATLAB 来求解非线性薛定谔方程的分步傅立叶数值算法。非线性薛定谔方程是由线性算符和非线性算符两部分组成，依据分步傅立叶数值算法的 思想，当光场每通过一段微小距离的时候，先发生线性算符代的色散数学方程，再计算非 线性算符代表的非线性效应数学方程，从而得到近似的数值结果。在计算存在高阶偏微分项 时，在时域中计算非常不便，可利用傅立叶变换，把偏微分方程变换为代数方程来运算，进 而求出其相应的数值结果。从算法的编写角度出发，该算法在 MATLAB 中易于实现，而且 比较简单清晰，能够广泛的应用在非线性薛定谔方程的数值求解之中。 关键词:非线性薛定谔方程;傅立叶变换;分步傅立叶数值算法;MATLAB 中图分类号:O 415 1. 引 言 非线性薛定谔方程是研究光脉冲在光纤中传输的基本方程，大多数文献都是直接引用非 线性薛定谔方程，采用分步傅里叶法数值求解非线性薛定谔方程。本文从通过MATLAB编 程来求解非线性薛定谔方程，并给出了分步傅里叶法数值求解算法程序。 MATLAB是美国MathWorks公司自20世纪80年代中期推出的数学软件，优秀的数值计算 能力和卓越的数据可视化能力使其很快在数学软件中脱颖而出。随着版本的不断升级，它在 数值计算及符号计算功能上得到了进一步完善。到目前为止，MATLAB已经发展成为多学 科、多种工作平台的功能强大的大型软件。 2. 非线性薛定谔方程 非线性薛定谔方程，简称 NLS 方程，是一个非线性偏微分方程，通常情况下无法求出 解析解，只能求出它的数值解。在含有非线性色散介质地脉冲传输问题中应用非常广泛。 非线性薛定谔方程的基本形式为: 2 (2-1) iu, u ， 2 | u | u txx 其中 u 是未知的复值函数. 3. 分步傅立叶数值算法 目前，采用分步傅立叶算法(Split step Fourier Method)求解非线性薛定谔方程的数值解应 用比较多。分步傅立叶方法最早是在1937年开始应用的，这种方法己经被证明是相同精度下 数值求解非线性薛定愕方程最快的方法，部分原因是它采用了快速傅立叶变换算法(Fast Fourier Transform Algorithm)。基于MATLAB科学计算软件以及MATLAB强大的符号计算功 能，完全可以实现分步傅立叶数值算法来对脉冲形状和频谱进行仿真。 - 1 - 如果您需要更多论文可以到www.docin.com/week114进行免费查阅。 一般情况下，光脉冲沿光纤传播时受到色散和非线性效应的共同作用，假设当传输距离很小的时候，两者相互独立作用，那么，根据这种思想可建立如下分步傅立叶数值算法的数 学模型: 把待求解的非线性薛定谔方程写成以下形式: ,U (3-1) ˆ ˆ , (D， N)U,z ˆ ˆ 其中 D是线性算符，代表介质的色散和损耗，N是非线性算符，它决定了脉冲传输过程中 光纤的非线性效应。 一般来讲，沿光纤的长度方向，色散和非线性是同时作用的。分步傅立叶法假设在传输 过程中，光场每通过一小段距离h ，色散和非线性效应可以分别作用，得到近似结果。也就 是说脉冲从 z 到 z ， h 的传输过程中分两步进行。第一步，只有非线性作用，方程(2)式中的 ˆ ˆ D, 0 ;第二步，再考虑线性作用，方程(2)式中的 N, 0 . 这样方程(2)在这两步中可分别简化为: ,U ˆ , D,U,z (3-2) ,U ˆ , N,U,z 得到了上面两个方程(3-2)，就可以分别求解非线性作用方程和线性作用方程，然后讨论分步傅立叶法的数值算法。 由于方程(3-2)是一个偏微分方程，需要通过傅立叶变换把偏 微分方程转换为代数方 程，进行运算。傅立叶变换的定义如下: ，,~ , F[U ( z,T )] , ( z,, ) , U U ( z,T ) exp(i,T )dT,,, , (3-3) , ，, ~~ 1,1 ,U ( z,, ) exp(,i,T )dT F [U ( z,, )] , U ( z, T ) ,, ,,, 2, 在计算 F[U ( z,T )] 时一般采用快速傅立叶变换(FFT)算。为了保证精度要求，一般还需要反复调整纵向传输步长 z 和横向脉冲取样点数 T 来保证计算精度。 4. 分步傅立叶数值算法的 MATLAB 实现 为了研究分步傅立叶法在 Matlab 上的实现，下面通过举例来说明 Matlab 在求解的非线 性薛定谔方程的优越性。 现待求解的非线性薛定谔方程如下: 2 ,i , A ,A ,2 ， A , , ,i | A |A , 0(4-1)2,z 24 ,T 其 中 ，是 光 场 慢 变 复 振 幅 ， z 是 脉 冲 沿 光 纤 传 播 的 距 离 ;A( z, T ) 是色散系数;, (1/ w , km) 是非线性系数; ( ps / km) T , t , ,z, ,, 1/ v , v 是群速度;1 1 g g , , (1 / km) ) 是 光 纤 损 耗 系 数 ， 它 与 用 分 贝 表 示 的 损 耗 系 数 , (dB / km) 的 关 系 为 : dB , , 4.343, . dB P，P是入射脉冲的峰值功率， 首先，可以将方程(4-1)归一化振幅:U , A( z,T ) / 0 0 此时方程(4-1)可改写为: - 2 - 如果您需要更多论文可以到www.docin.com/week114进行免费查阅。 2 ,U , ,i , U2 , , U ， ， ,P i | U |U(4-2)02 ,z 24 ,T 为了使用分步傅立叶法求解方程(4-2)，将方程(4-2)写成以下形式: ,U ˆ ˆ , (D， N)U ,z 进一步，可以得出如下方程(4-3): 2,i , U , , ，2 2 ,T ˆ (4-3) D, 2 2ˆ N, ,P i | U |0 然后，按照步骤 step1 和步骤 step2，依次计算方程(4-3)的线性算符和非线性算符。最后在步骤 step3 中，运行步骤 step1 和步骤 step2 的 MATLAB 程序，得出线性算符和非线性算 符的精确数值解及其仿真曲线。 Step1 线性算符方程的求解 ˆ 线性算符 D的方程如下: 2,i , U , , ，2 ,U 2 ,T (4-4), U,z 2 用傅立叶变换方法，得到一个常微分方程(4-5): ~2 , , i(i, ) U ,~ ~ (4-5), , U , U 4 ,z 2 解方程(4-5)得: 2 ,, , ,~~iU U ( z, , ) ,(0,, ) exp[ (4-6) 2 z] 4 ~~ 式中 (0,, ) 是初值U (0, T ) 的傅立叶变换，将 ( z, ,) 进行反傅立叶变换就得到了U U U ( z, T ) 。方程(4-6)的求解公式为: i~ z 2 (4-7) F{e U ( z, T ) , xp[( ,, , , ) ] , F[U (0,T )]} 22 ~其中 F 和 F 分别表示傅立叶变换和反傅立叶变换运算。 Step2 非线性算符方程的求解 ˆ 非线性部分 N的方程如下: ,U 2 (4-8) , ,P i | U |U 0 ,z 同Step1的方法，解方程(4-8)，得到: ~~2 U U ( z, , ) ,(0,, ) exp[,P i | U (0, T ) |z](4-9) 0 ~~式中 (0,, ) 是初值U (0, T ) 的傅立叶变换，将 ( z, ,) 进行反傅立叶变换就得到了U U U ( z, T ) 。方程(4-9)的求解公式为: ~2 U ( z, T ) , F{exp[P ,i | U (0,T ) | z] , F[U (0,T )]} (4-10)0 - 3 - 如果您需要更多论文可以到www.docin.com/week114进行免费查阅。 ~ F 其中 F 和 分别表示傅立叶变换和反傅立叶变换运算。 Step 3 算法在MATLAB中的实现 在Matlab中，设有限时长序列 x(n) 的长度为 N (1 , n , N ) ，它对应于一个频域内的长 2,k (1 , k , N ) ，其中T 度为N的有限长序列 X (k )(1 , k , N ) . X (k ) 对应的角频, (k ) , NT 是序列 x(n) 的采样时间间隔. 这种正反离散傅立叶变换的关系式为: N 2, X (k ) , DFT [ x(n)] , x(n) exp(, j , , k , n)(1 , k , N ), N j ,1 (4-11) N1 2, x(n) , IDFT [ X (k )] ,X (k ) exp( j , , k , n)(1 , n , N ), N N j ,1 然后用Matlab中的离散傅立叶变换(DFT)函数fft和离散傅立叶反变换(IDFT)的函数ifft来 实现方程(4-4)，(4-8)式中的傅立叶和反傅立叶变换运算。进一步，得到方程(4-7)，(4-10)的 数值解及仿真曲线。 实现分步傅立叶变换算法MATLAB的脚本文件程序如下: Po;%输入光强，单位W alpha;%光纤损耗值，单位dB/km gamma;%光纤非线性参数 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?%步长 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))? - 4 - 如果您需要更多论文可以到www.docin.com/week114进行免费查阅。 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:z spectrum=spectrum.*exp(g1) ? %g1为线性算符e的指数表达式 m)? %快速离散反傅立叶变换 f=ifft(spectru f=f.*exp(g2)?%g2为非线性算符e的指数表达式 spectrum=fft(f)? %快速离散傅立叶变换 spectrum=spectrum.*exp(g1) ? end 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? =atand((abs(imag(f)))/(abs(real(f))))? dd phadisp(ln)=dd?%保存脉冲相位 ln=ln+1? end toc? 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(5) plot(phadisp(1:1:ln-1),'k')? xlabel('distance travelled')? ylabel('phase change')? grid on? hold on? disp('CPU time:'), disp(cputime)? 具体在实现分步傅立叶变换算法的时候，除了上述脚本本件定义的参数外，还可以根据 实际情况，修改程序部分参数，直到达到你想要求解精度和曲线。 最后，通过测试一组参 数，得到方程(4-1)在该算法下的MATLAB运算结果。MATLAB 总共用时177.2660s，求得的的结果曲线如下图1所示。结果表明，算法正确而且精度也比较 高，能够在非线性薛定谔方程的求解中广泛应用。 - 5 - 如果您需要更多论文可以到www.docin.com/week114进行免费查阅。 图1. 仿真结果 参考文献 [1] 程雪苹. 直接微扰方法和非线性薛定谔方程族的近似解[D]. 浙江:浙江师范大学,2006. [2] 齐志勇. 掺铒光纤放大器的模拟器研究[D]. 湖北:华中科技大学,2004. [3] 施娟. 基于对称分步傅立叶算法的光孤子仿真,J,. 技术前沿,2008,第 10 卷(第 1 期) [4] 张兴坊,郑义,李爱萍,徐云峰.皮秒脉冲在光子晶体光纤中的压缩效应,J,. 激光技术,200,第 31 卷(第 3 期) [5] 吕理想,张晓萍.不同形式非线性薛定谔方程及其分步傅里叶法求解,J,. 计算物理,2007,第 24 卷(第 3 期) [6] 张志涌.《精通 MATLAB6.5 版》,M,. 北京:北京航空航天大学出版社,2003.. Solutions for nonlinear Schrödinger equation by Split step Fourier Method on MATLAB LEI Chao Colledge of Agricultural Engeering, Sichuan Agricultural University Abstract - 6 - 如果您需要更多论文可以到www.docin.com/week114进行免费查阅。 Schrödinger Equation can be solved by Split Step Fourier Method paper,Nonlinear In this Algorithm based on MATLAB. Nonlinear Schrödinger Equation is composed of linear functor and nonlinear functor. According to the theory of Split Step Fourier Method Algorithm, the linear mathematical equation will be first worked out,and then calculate the nonlinear mathematical equation when the LightField pass a slight distance. If there exists the high-level order partial differential equation which cannot be calculated in time-domain, by the method of Fourier Transform, the high- level order partial differential equation will be transformed to the algebra equation. The numerical value results will be concluded by this way. From compiling the algorithm point of view, this algorithm is simple,clear and easy to be carried out on MATLAB, and that is able to apply extensively for the solutions of nonlinear Schrödinger equation. Keywords: Nonlinear Schrödinger Equation, Fourier Transform, Split Step Fourier Method Algorithm, MATLAB - 7 -