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

首先， 可以将方程 （4­1） 归一化振幅：U A( z , T ) /
P0 ，P0 是入射脉冲的峰值功率，
此时方程(4­1)可改写为:
­ 2 ­

U i 2U U P0 i | U | 2 U 2 z 2 4 T
（4­2）
3. 分步傅立叶数值算法
目前， 采用分步傅立叶算法(Split step Fourier Method)求解非线性薛定谔方程的数值解应 用比较多。 分步傅立叶方法最早是在1937年开始应用的， 这种方法己经被证明是相同精度下 数值求解非线性薛定愕方程最快的方法，部分原因是它采用了快速傅立叶变换算法 (Fast Fourier Transform Algorithm)。基于MATLAB科学计算软件以及MATLAB强大的符号计算功 能，完全可以实现分步傅立叶数值算法来对脉冲形状和频谱进行仿真。
其中 F 和 F 分别表示傅立叶变换和反傅立叶变换运算。 Step2 非线性算符方程的求解
~
i z ~ U ( z , T ) F{exp[( 2 ) ] F [U (0, T )]} 2 2
Hale Waihona Puke （4­7）ˆ 的方程如下： 非线性部分 N
U P0 i | U | 2 U z
为了使用分步傅立叶法求解方程(4­2)，将方程(4­2)写成以下形式:
U ˆN ˆ )U (D z
进一步，可以得出如下方程（4­3）：
i 2U 2 T 2 ˆ D 2 ˆ P i | U | 2 N 0
（4­3）
然后，按照步骤 step1 和步骤 step2，依次计算方程（4­3）的线性算符和非线性算符。最后 在步骤 step3 中，运行步骤 step1 和步骤 step2 的 MATLAB 程序，得出线性算符和非线性算 符的精确数值解及其仿真曲线。 Step1 线性算符方程的求解
ˆ 0 ；第二步，再考虑线性作用，方程(2)式中的 N ˆ 0. D
这样方程(2)在这两步中可分别简化为:
U ˆ U D z U ˆ U N z
分步傅立叶法的数值算法。
（3­2）
得到了上面两个方程（3­2），就可以分别求解非线性作用方程和线性作用方程，然后讨论 由于方程（3­2）是一个偏微分方程，需要通过傅立叶变换把偏微分方程转换为代数方 程，进行运算。傅立叶变换的定义如下：
2. 非线性薛定谔方程
非线性薛定谔方程，简称 NLS 方程，是一个非线性偏微分方程，通常情况下无法求出 解析解，只能求出它的数值解。在含有非线性色散介质地脉冲传输问题中应用非常广泛。 非线性薛定谔方程的基本形式为:
iu t u xx 2 | u | 2 u
其中 u 是未知的复值函数.
（2­1）
光纤的非线性效应。
（3­1）
ˆ 是线性算符，代表介质的色散和损耗， N ˆ 是非线性算符，它决定了脉冲传输过程中 其中 D
一般来讲，沿光纤的长度方向，色散和非线性是同时作用的。分步傅立叶法假设在传输 过程中，光场每通过一小段距离 h ，色散和非线性效应可以分别作用，得到近似结果。也就 是说脉冲从 z 到 z h 的传输过程中分两步进行。第一步，只有非线性作用，方程(2)式中的
1. 引 言
非线性薛定谔方程是研究光脉冲在光纤中传输的基本方程， 大多数文献都是直接引用 非 线性薛定谔方程，采用分步傅里叶法数值求解非线性薛定谔方程。本文从通过MATLAB编 程来求解非线性薛定谔方程，并给出了分步傅里叶法数值求解算法程序。 MATLAB是美国MathWorks公司自20世纪80年代中期推出的数学软件， 优秀的数值计算 能力和卓越的数据可视化能力使其很快在数学软件中脱颖而出。 随着版本的不断升级， 它 在 数值计算及符号计算功能上得到了进一步完善。到目前为止，MATLAB已经发展成为多学 科、多种工作平台的功能强大的大型软件。
反复调整纵向传输步长 z 和横向脉冲取样点数 T 来保证计算精度。
4. 分步傅立叶数值算法的 MATLAB 实现
为了研究分步傅立叶法在 Matlab 上的实现，下面通过举例来说明 Matlab 在求解的非线 性薛定谔方程的优越性。 现待求解的非线性薛定谔方程如下:
A i 2 A A i | A | 2 A 0 z 2 4 T 2
（4­8）
同Step1的方法，解方程（4­8），得到： ~ ~ U ( z , ) U (0, ) exp[P0 i | U (0, T ) | 2 z ] （4­9） ~ ~ 式中 U (0, ) 是初值 U (0, T ) 的傅立叶变换，将U ( z , ) 进行反傅立叶变换就得到了 U ( z , T ) 。方程(4­9)的求解公式为: ~ U ( z , T ) F{exp[ P0i | U (0, T ) | 2 z ] F [U (0, T )]} （4­10）
~
2k (1 k N ) ，其中 T NT
是序列 x(n) 的采样时间间隔. 这种正反离散傅立叶变换的关系式为:
N
X (k ) DFT [ x(n)] x(n) exp( j
j 1
2 k n)(1 k N ) N
1 x(n) IDFT [ X (k )] N
­ 3 ­

其中 F 和 F 分别表示傅立叶变换和反傅立叶变换运算。 Step 3 算法在MATLAB中的实现 在Matlab中，设有限时长序列 x(n) 的长度为 N (1 n N ) ，它对应于一个频域内的长 度为N的有限长序列 X ( k )(1 k N ) . X ( k ) 对应的角频 ( k )
ˆ 的方程如下： 线性算符 D
U z

i 2U 2 T 2 U 2
（4­4）
用傅立叶变换方法，得到一个常微分方程（4­5）：
~ U ~ i (i ) 2 ~ U U z 2 4
解方程(4­5)得:
（4­5）
i 2 2 ~ ~ （4­6） U ( z , ) U (0, ) exp[ z] 4 ~ ~ 式中 U (0, ) 是初值 U (0, T ) 的傅立叶变换，将U ( z , ) 进行反傅立叶变换就得到了 U ( z , T ) 。方程(4­6)的求解公式为:
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)); %扩散长度，单位是m Ao=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));
­ 1 ­

一般情况下，光脉冲沿光纤传播时受到色散和非线性效应的共同作用，假设当传输距离 很小的时候，两者相互独立作用，那么，根据这种思想可建立如下分步傅立叶数值算法的数 学模型： 把待求解的非线性薛定谔方程写成以下形式:
U ˆN ˆ )U (D z
（4­1）
其 中 ， A( z , T ) 是 光 场 慢 变 复 振 幅 ， z 是 脉 冲 沿 光 纤 传 播 的 距 离 ；
(1 / w km) 是非线性系数； T t 1 z , 1 1 / v g , v g 是群速度； ( ps / km) 是色散系数； (1 / km) ) 是 光 纤 损 耗 系 数 ， 它 与 用 分 贝 表 示 的 损 耗 系 数 dB (dB / km) 的 关 系 为 ： dB 4.343 .
2 X (k ) exp( j k n)(1 n N ) N j 1
N
（4­11）
然后用Matlab中的离散傅立叶变换(DFT)函数fft和离散傅立叶反变换(IDFT)的函数ifft来 实现方程(4­4)，(4­8)式中的傅立叶和反傅立叶变换运算。进一步，得到方程(4­7)，(4­10)的 数值解及仿真曲线。 实现分步傅立叶变换算法MATLAB的脚本文件程序如下：
~ F [U ( z , T )] U ( z , ) U ( z , T ) exp(iT )dT （3­3） 1 ~ 1 ~ F [U ( z , )] U ( z , T ) U ( z , ) exp(iT )dT 2 在计算 F [U ( z , T )] 时一般采用快速傅立叶变换（FFT）算。为了保证精度要求，一般还需要
­ 4 ­
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的指数表达式 f=ifft(spectrum); %快速离散反傅立叶变换 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; dd=atand((abs(imag(f)))/(abs(real(f)))); 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);