非線性薛定諤方程數(shù)值解的MATLAB仿真

1、admin鍵入文檔標題利用分步快速傅里葉變換對光纖中光信號的傳輸方程進行數(shù)值求解1、非線性薛定諤方程非線性薛定諤方程(nonlinear Schrodinger equation，NLSE)是奧地利物理學家薛定諤于1926 年提出的，應用在量子力學系統(tǒng)中。由于量子力學主要研究粒子的動力學運動狀態(tài)，所以不能運用牛頓力學公式來表示。通常在量子力學中，研究系統(tǒng)的狀態(tài)一般通過波函數(shù)(x，t)來表示。而對波函數(shù)的研究主要是求解非線性薛定諤方程。本文主要研究光脈沖在光纖中傳輸狀態(tài)下的演變。一般情況下，光脈沖信號在光纖中傳輸時，同時受到光纖的色散和非線性效應的影響。通過Maxwell 方程，考慮到光纖的色散

2、和非線性效應，可以推導出光信號在光纖中的傳輸方程，即非線性薛定諤方程。NLSE 是非線性偏微分方程，一般很難直接求出解析解，于是通過數(shù)值方法進行求解。具體分為兩大類：(1)分布有限差分法(split-step finite differencemethod，SSFD)；(2)分步傅里葉變換法(split-step Fourier transform method，SSFT)。一般情況，在達到相同精度，由于分步傅里葉變換法采用運算速度快的快速傅里葉變換，所以相比較有限差分法運算速度快一到兩個數(shù)量級。于是本文介紹分步傅里葉變換法來對光纖中光信號的傳輸方程，即非線性薛定諤方程進行數(shù)值求解。并通過MA

3、TLAB 軟件對結果數(shù)值仿真。非線性薛定諤方程的基本形式為:（I）其中u 是未知的復值函數(shù).目前，采用分步傅立葉算法(Split step Fourier Method)求解非線性薛定諤方程的數(shù)值解應用比較多。分步傅立葉方法最早是在1937年開始應用的，這種方法己經被證明是相同精度下數(shù)值求解非線性薛定愕方程最快的方法，部分原因是它采用了快速傅立葉變換算法(Fast Fourier Transform Algorithm)。基于MATLAB科學計算軟件以及MATLAB強大的符號計算功能，完全可以實現(xiàn)分步傅立葉數(shù)值算法來對脈沖形狀和頻譜進行仿真。一般情況下，光脈沖沿光纖傳播時受到色散和非線性效應的

4、共同作用，假設當傳輸距離很小的時候，兩者相互獨立作用，那么，根據(jù)這種思想可建立如下分步傅立葉數(shù)值算法的數(shù)學模型：把待求解的非線性薛定諤方程寫成以下形式:（II）¶其中是線性算符，代表介質的色散和損耗， 是非線性算符，它決定了脈沖傳輸過程中光纖的非線性效應。一般來講，沿光纖的長度方向，色散和非線性是同時作用的。分步傅立葉法假設在傳輸過程中，光場每通過一小段距離h，色散和非線性效應可以分別作用，得到近似結果。也就是說脈沖從z到z+h的傳輸過程中分兩步進行。第一步，只有非線性作用，方程(II)式中的=0；第二步，再考慮線性作用，方程(II)式中的=0這樣方程(2)在這兩步中可分別簡化為:（

5、III）得到了上面兩個方程（III），就可以分別求解非線性作用方程和線性作用方程，然后討論分步傅立葉法的數(shù)值算法。由于方程（III）是一個偏微分方程，需要通過傅立葉變換把偏微分方程轉換為代數(shù)方程，進行運算。傅立葉變換的定義如下：（IV）在計算時一般采用快速傅立葉變換（FFT）算。為了保證精度要求，一般還需要反復調整縱向傳輸步長z和橫向脈沖取樣點數(shù)T來保證計算精度。2、分步傅立葉數(shù)值算法的MATLAB 實現(xiàn)現(xiàn)待求解的非線性薛定諤方程如下:（V）其中，A(z,T)是光場慢變復振幅，z是脈沖沿光纖傳播的距離；,vg是群速度；是色散系數(shù)；是非線性系數(shù)；是光纖損耗系數(shù)，它與用分貝表示的損耗系數(shù)的關系為：

6、.首先，可以將方程（V）歸一化振幅：， 是入射脈沖的峰值功率，此時方程(V)可改寫為：（VI）為了使用分步傅立葉法求解方程(VI)，將方程(VI)寫成以下形式:進一步，可以得出如下方程（VII）：（VII）然后，按照步驟1和步驟2，依次計算方程（VII）的線性算符和非線性算符。最后在步驟3 中，運行步驟1和步驟2的MATLAB 程序，得出線性算符和非線性算符的精確數(shù)值解及其仿真曲線。步驟1 線性算符方程的求解線性算符的方程如下： （VIII） ba用傅立葉變換方法，得到一個常微分方程（IX）：（IX）解方程(IX)得:（X）式中是初值的傅立葉變換，將進行反傅立葉變換就得到了。方程(X)的求解公

7、式為:（XI）其中和分別表示傅立葉變換和反傅立葉變換運算。步驟2 非線性算符方程的求解非線性部分的方程如下： （XII）同Step1的方法，解方程（XII），得到：（XIII）式中是初值的傅立葉變換，將進行反傅立葉變換就得到了。方程(XIII)的求解公式為:（XIV）其中和分別表示傅立葉變換和反傅立葉變換運算。步驟3 算法在MATLAB中的實現(xiàn)在Matlab中，設有限時長序列的長度為，它對應于一個頻域內的長度為N的有限長序列，的角頻，其中T是序列的采樣時間間隔.這種正反離散傅立葉變換的關系式為:（XV）然后用Matlab中的離散傅立葉變換(DFT)函數(shù)fft和離散傅立葉反變換(IDFT)的函數(shù)

8、ifft來實現(xiàn)方程(VIII)，(XII)式中的傅立葉和反傅立葉變換運算。進一步，得到方程(XI)，(XIV)的數(shù)值解及仿真曲線。最后，通過測試一組參數(shù)，得到方程（V）在該算法下的MATLAB運算結果。MATLAB總共用時34.26s，求得的的結果曲線如下圖所示。結果表明，算法正確而且精度也比較高，能夠在非線性薛定諤方程的求解中廣泛應用。附錄MATLAB的腳本文件源代碼：Po=200; %輸入光強，單位Walpha=3.5; %光纖損耗值，單位dB/kmgamma=20; %光纖非線性參數(shù)to=1; %初始脈沖寬度，單位秒C=50; %第一次計算輸入的啁啾參數(shù)b2=1000; %波數(shù)的倒數(shù)cp

9、utime=0;tic;ln=1;i=sqrt(-1);pi=3.1415926535;alph=alpha/(4.343);Ld=(to2)/(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'

10、); 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=(-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); %

11、g1為線性算符e的指數(shù)表達式f=ifft(spectrum); %快速離散反傅立葉變換f=f.*exp(g2);%g2為非線性算符e的指數(shù)表達式spectrum=fft(f); %快速離散傅立葉變換spectrum=spectrum.*exp(g1);endf=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(

12、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 of steps');ylabel('Pulse broadening ratio&#