非线性薛定谔方程数值解的MATLAB仿真_精品文档Word文档下载推荐.docx

非线性薛定谔方程数值解的MATLAB仿真_精品文档Word文档下载推荐.docx

《非线性薛定谔方程数值解的MATLAB仿真_精品文档Word文档下载推荐.docx》由会员分享，可在线阅读，更多相关《非线性薛定谔方程数值解的MATLAB仿真_精品文档Word文档下载推荐.docx（9页珍藏版）》请在冰豆网上搜索。

其中u是未知的复值函数.

目前，采用分步傅立叶算法（SplitstepFourierMethod）求解非线性薛定谔方程的数值解应用比较多。

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

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

一般情况下，光脉冲沿光纤传播时受到色散和非线性效应的共同作用，假设当传输距离

很小的时候，两者相互独立作用，那么，根据这种思想可建立如下分步傅立叶数值算法的数

学模型：

把待求解的非线性薛定谔方程写成以下形式:

（II）

¶

其中是线性算符，代表介质的色散和损耗，是非线性算符，它决定了脉冲传输过程中光纤的非线性效应。

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

分步傅立叶法假设在传输

过程中，光场每通过一小段距离h，色散和非线性效应可以分别作用，得到近似结果。

也就

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

第一步，只有非线性作用，方程（II）式中的

=0；

第二步，再考虑线性作用，方程（II）式中的=0

这样方程

（2）在这两步中可分别简化为:

（III）

得到了上面两个方程（III），就可以分别求解非线性作用方程和线性作用方程，然后讨论分步傅立叶法的数值算法。

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

程，进行运算。

傅立叶变换的定义如下：

（IV）

在计算时一般采用快速傅立叶变换（FFT）算。

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

2、分步傅立叶数值算法的MATLAB实现

现待求解的非线性薛定谔方程如下:

（V）

其中，A（z,T）是光场慢变复振幅，z是脉冲沿光纤传播的距离；

,vg是群速度；

是色散系数；

是非线性系数；

是光纤损耗系数，它与用分贝表示的损耗系数的关系为：

.

首先，可以将方程（V）归一化振幅：

，是入射脉冲的峰值功率，

此时方程（V）可改写为：

（VI）

为了使用分步傅立叶法求解方程（VI），将方程（VI）写成以下形式:

进一步，可以得出如下方程（VII）：

（VII）

然后，按照步骤1和步骤2，依次计算方程（VII）的线性算符和非线性算符。

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

步骤1线性算符方程的求解

线性算符的方程如下：

（VIII）

b

a

用傅立叶变换方法，得到一个常微分方程（IX）：

（IX）

解方程（IX）得:

（X）

式中是初值的傅立叶变换，将进行反傅立叶变换就得到了。

方程（X）的求解公式为:

（XI）

其中和分别表示傅立叶变换和反傅立叶变换运算。

步骤2非线性算符方程的求解

非线性部分的方程如下：

（XII）

同Step1的方法，解方程（XII），得到：

（XIII）

方程（XIII）的求解公式为:

（XIV）

步骤3算法在MATLAB中的实现

在Matlab中，设有限时长序列的长度为，它对应于一个频域内的长度为N的有限长序列，的角频，其中T是序列的采样时间间隔.

这种正反离散傅立叶变换的关系式为:

（XV）

然后用Matlab中的离散傅立叶变换（DFT）函数fft和离散傅立叶反变换（IDFT）的函数ifft来实现方程（VIII），（XII）式中的傅立叶和反傅立叶变换运算。

进一步，得到方程（XI），（XIV）

的数值解及仿真曲线。

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

MATLAB总共用时34.26s，求得的的结果曲线如下图所示。

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

附录

MATLAB的脚本文件源代码：

Po=200;

%输入光强，单位W

alpha=3.5;

%光纤损耗值，单位dB/km

gamma=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））;

%扩散长度，单位是m

Ao=sqrt（Po）;

%光振幅

tau=-4096e-12:

1e-12:

4095e-12;

dt=1e-12;

h=1000;

%步长

forii=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（'

InputPulse'

xlabel（'

Time'

ylabel（'

Amplitude'

gridon;

holdon;

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））;

%快速离散傅立叶变换

forjj=h:

h:

z

spectrum=spectrum.*exp（g1）;

%g1为线性算符e的指数表达式

f=ifft（spectrum）;

%快速离散反傅立叶变换

f=f.*exp（g2）;

%g2为非线性算符e的指数表达式

spectrum=fft（f）;

end

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;

toc;

cputime=toc;

figure

（2）;

mesh（op_pulse（1:

ln-1,:

））;

title（'

PulseEvolution'

xlabel（'

distance'

zlabel（'

amplitude'

figure（3）

plot（pbratio（1:

ln-1）,'

k'

Numberofsteps'

ylabel（'

Pulsebroadeningratio'

figure（4）

plot（phadisp（1:

distancetravelled'

phasechange'

disp（'

CPUtime:

'

）,disp（cputime）;