基于matlab的非线性薛定谔方程的数值算法研究
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
首先， 可以将方程 （41） 归一化振幅：U A( z , T ) /
P0 ，P0 是入射脉冲的峰值功率，
此时方程(41)可改写为:
2
U i 2U U P0 i | U | 2 U 2 z 2 4 T
（42）
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 （47）ˆ 的方程如下： 非线性部分 N
U P0 i | U | 2 U z
为了使用分步傅立叶法求解方程(42)，将方程(42)写成以下形式:
U ˆN ˆ )U (D z
进一步，可以得出如下方程（43）：
i 2U 2 T 2 ˆ D 2 ˆ P i | U | 2 N 0
（43）
然后，按照步骤 step1 和步骤 step2，依次计算方程（43）的线性算符和非线性算符。最后 在步骤 step3 中，运行步骤 step1 和步骤 step2 的 MATLAB 程序，得出线性算符和非线性算 符的精确数值解及其仿真曲线。 Step1 线性算符方程的求解
ˆ 0 ；第二步，再考虑线性作用，方程(2)式中的 N ˆ 0. D
这样方程(2)在这两步中可分别简化为:
U ˆ U D z U ˆ U N z
分步傅立叶法的数值算法。
（32）
得到了上面两个方程（32），就可以分别求解非线性作用方程和线性作用方程，然后讨论 由于方程（32）是一个偏微分方程，需要通过傅立叶变换把偏微分方程转换为代数方 程，进行运算。傅立叶变换的定义如下：
2. 非线性薛定谔方程
非线性薛定谔方程，简称 NLS 方程，是一个非线性偏微分方程，通常情况下无法求出 解析解，只能求出它的数值解。在含有非线性色散介质地脉冲传输问题中应用非常广泛。 非线性薛定谔方程的基本形式为:
iu t u xx 2 | u | 2 u
其中 u 是未知的复值函数.
（21）
光纤的非线性效应。
（31）
ˆ 是线性算符，代表介质的色散和损耗， 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
（48）
同Step1的方法，解方程（48），得到： ~ ~ U ( z , ) U (0, ) exp[P0 i | U (0, T ) | 2 z ] （49） ~ ~ 式中 U (0, ) 是初值 U (0, T ) 的傅立叶变换，将U ( z , ) 进行反傅立叶变换就得到了 U ( z , T ) 。方程(49)的求解公式为: ~ U ( z , T ) F{exp[ P0i | U (0, T ) | 2 z ] F [U (0, T )]} （410）
~
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
（44）
用傅立叶变换方法，得到一个常微分方程（45）：
~ U ~ i (i ) 2 ~ U U z 2 4
解方程(45)得:
（45）
i 2 2 ~ ~ （46） U ( z , ) U (0, ) exp[ z] 4 ~ ~ 式中 U (0, ) 是初值 U (0, T ) 的傅立叶变换，将U ( z , ) 进行反傅立叶变换就得到了 U ( z , T ) 。方程(46)的求解公式为:
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 = 4096e12:1e12: 4095e12; dt=1e12; 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
（41）
其 中 ， 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
（411）
然后用Matlab中的离散傅立叶变换(DFT)函数fft和离散傅立叶反变换(IDFT)的函数ifft来 实现方程(44)，(48)式中的傅立叶和反傅立叶变换运算。进一步，得到方程(47)，(410)的 数值解及仿真曲线。 实现分步傅立叶变换算法MATLAB的脚本文件程序如下：
~ F [U ( z , T )] U ( z , ) U ( z , T ) exp(iT )dT （33） 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/21)*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:ln1,:)); title('Pulse Evolution'); xlabel('Time'); ylabel('distance'); zlabel('amplitude'); figure(3) plot(pbratio(1:1:ln1),'k'); xlabel('Number of steps'); ylabel('Pulse broadening ratio'); grid on; hold on; figure(5) plot(phadisp(1:1:ln1),'k'); xlabel('distance travelled'); ylabel('phase change'); grid on; hold on; disp('CPU time:'), disp(cputime);
基于 MATLAB 的非线性薛定谔方程的数值算法研究
雷超
四川农业大学农业工程系，四川雅安 （625014） Email（leichaocool123@）
摘 要：本文简单介绍了一种基于 MATLAB 来求解非线性薛定谔方程的分步傅立叶数值算 法。 非线性薛定谔方程是由线性算符和非线性算符两部分组成， 依据分步傅立叶数值算法的 思想，当光场每通过一段微小距离的时候，先发生线性算符代表的色散数学方程，再计算非 线性算符代表的非线性效应数学方程， 从而得到近似的数值结果。 在计算存在高阶偏微分项 时，在时域中计算非常不便，可利用傅立叶变换，把偏微分方程变换为代数方程来运算， 进 而求出其相应的数值结果。从算法的编写角度出发，该算法在 MATLAB 中易于实现，而且 比较简单清晰，能够广泛的应用在非线性薛定谔方程的数值求解之中。 关键词：非线性薛定谔方程；傅立叶变换；分步傅立叶数值算法；MATLAB 中图分类号：O 415
P0 ，P0 是入射脉冲的峰值功率，
此时方程(41)可改写为:
2
U i 2U U P0 i | U | 2 U 2 z 2 4 T
（42）
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 （47）ˆ 的方程如下： 非线性部分 N
U P0 i | U | 2 U z
为了使用分步傅立叶法求解方程(42)，将方程(42)写成以下形式:
U ˆN ˆ )U (D z
进一步，可以得出如下方程（43）：
i 2U 2 T 2 ˆ D 2 ˆ P i | U | 2 N 0
（43）
然后，按照步骤 step1 和步骤 step2，依次计算方程（43）的线性算符和非线性算符。最后 在步骤 step3 中，运行步骤 step1 和步骤 step2 的 MATLAB 程序，得出线性算符和非线性算 符的精确数值解及其仿真曲线。 Step1 线性算符方程的求解
ˆ 0 ；第二步，再考虑线性作用，方程(2)式中的 N ˆ 0. D
这样方程(2)在这两步中可分别简化为:
U ˆ U D z U ˆ U N z
分步傅立叶法的数值算法。
（32）
得到了上面两个方程（32），就可以分别求解非线性作用方程和线性作用方程，然后讨论 由于方程（32）是一个偏微分方程，需要通过傅立叶变换把偏微分方程转换为代数方 程，进行运算。傅立叶变换的定义如下：
2. 非线性薛定谔方程
非线性薛定谔方程，简称 NLS 方程，是一个非线性偏微分方程，通常情况下无法求出 解析解，只能求出它的数值解。在含有非线性色散介质地脉冲传输问题中应用非常广泛。 非线性薛定谔方程的基本形式为:
iu t u xx 2 | u | 2 u
其中 u 是未知的复值函数.
（21）
光纤的非线性效应。
（31）
ˆ 是线性算符，代表介质的色散和损耗， 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
（48）
同Step1的方法，解方程（48），得到： ~ ~ U ( z , ) U (0, ) exp[P0 i | U (0, T ) | 2 z ] （49） ~ ~ 式中 U (0, ) 是初值 U (0, T ) 的傅立叶变换，将U ( z , ) 进行反傅立叶变换就得到了 U ( z , T ) 。方程(49)的求解公式为: ~ U ( z , T ) F{exp[ P0i | U (0, T ) | 2 z ] F [U (0, T )]} （410）
~
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
（44）
用傅立叶变换方法，得到一个常微分方程（45）：
~ U ~ i (i ) 2 ~ U U z 2 4
解方程(45)得:
（45）
i 2 2 ~ ~ （46） U ( z , ) U (0, ) exp[ z] 4 ~ ~ 式中 U (0, ) 是初值 U (0, T ) 的傅立叶变换，将U ( z , ) 进行反傅立叶变换就得到了 U ( z , T ) 。方程(46)的求解公式为:
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 = 4096e12:1e12: 4095e12; dt=1e12; 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
（41）
其 中 ， 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
（411）
然后用Matlab中的离散傅立叶变换(DFT)函数fft和离散傅立叶反变换(IDFT)的函数ifft来 实现方程(44)，(48)式中的傅立叶和反傅立叶变换运算。进一步，得到方程(47)，(410)的 数值解及仿真曲线。 实现分步傅立叶变换算法MATLAB的脚本文件程序如下：
~ F [U ( z , T )] U ( z , ) U ( z , T ) exp(iT )dT （33） 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/21)*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:ln1,:)); title('Pulse Evolution'); xlabel('Time'); ylabel('distance'); zlabel('amplitude'); figure(3) plot(pbratio(1:1:ln1),'k'); xlabel('Number of steps'); ylabel('Pulse broadening ratio'); grid on; hold on; figure(5) plot(phadisp(1:1:ln1),'k'); xlabel('distance travelled'); ylabel('phase change'); grid on; hold on; disp('CPU time:'), disp(cputime);
基于 MATLAB 的非线性薛定谔方程的数值算法研究
雷超
四川农业大学农业工程系，四川雅安 （625014） Email（leichaocool123@）
摘 要：本文简单介绍了一种基于 MATLAB 来求解非线性薛定谔方程的分步傅立叶数值算 法。 非线性薛定谔方程是由线性算符和非线性算符两部分组成， 依据分步傅立叶数值算法的 思想，当光场每通过一段微小距离的时候，先发生线性算符代表的色散数学方程，再计算非 线性算符代表的非线性效应数学方程， 从而得到近似的数值结果。 在计算存在高阶偏微分项 时，在时域中计算非常不便，可利用傅立叶变换，把偏微分方程变换为代数方程来运算， 进 而求出其相应的数值结果。从算法的编写角度出发，该算法在 MATLAB 中易于实现，而且 比较简单清晰，能够广泛的应用在非线性薛定谔方程的数值求解之中。 关键词：非线性薛定谔方程；傅立叶变换；分步傅立叶数值算法；MATLAB 中图分类号：O 415