計算物理實習(xí)

1、 計 算 物 理 實 習(xí)即應(yīng)用 MATLAB仿真研究孤子對在飽和介質(zhì)中的傳輸 （0640502117 余方玉） （0640502110 李哲軒）研究文章： 孤子對在飽和介質(zhì)中的傳輸文章作者： 方云團 王永順 沈廷根 李利方摘要：孤子對在飽和介質(zhì)中的傳輸一文中，作者利用分步傅里葉變換法求解了孤子對在飽和非線性介質(zhì)中傳輸?shù)腘LS方程，本文在此解的基礎(chǔ)上利用MATLAB仿真得到在此條件下孤子對傳輸?shù)臄?shù)值圖形。關(guān)鍵詞：孤子對傳輸；MATLAB仿真1、仿真公式來源非線性薛定諤方程（NLS）是研究光孤子通信的動力學(xué)基礎(chǔ)，但該方程所規(guī)定的非線性大多屬于Kerr類型，而當(dāng)注入光脈沖能量較大時，必須考慮介質(zhì)折射

2、率的非線性飽和效應(yīng)，無耗光纖中孤子的場包絡(luò)函數(shù)q(x,)的傳輸方程為： （1）P為色散系數(shù)，且p>0,A為光纖非線性系數(shù)，H為描述介質(zhì)非線性飽和因子，文中作者在計算時取p=0.5,A=1.分步傅里葉變換后得： （2）上面公式中的和就是傅里葉變換和逆變換。2、仿真步驟和結(jié)果 我們?nèi)〕跏济}沖 ，其圖形如下 圖1：初始脈沖 MATLAB中提供了fft和ifft函數(shù)進行上述傅里葉變換及其逆變換，在數(shù)值計算時我們也是主要應(yīng)用這兩個函數(shù)，在公式（2）中我們看到在傅里葉變換之后要乘以，里面的是頻域空間的變量，所以計算時需要在傅里葉變換之后把時域空間變到頻域空間，在此我們應(yīng)用wspace函數(shù)來完成這個變

3、換。公式（2）中的H=0時數(shù)值解結(jié)果如圖2和圖3： 圖2：H=0時孤子對在飽和介質(zhì)中的傳輸 圖3：H=0時孤子對傳輸?shù)慕孛鎴D： 公式（2）中H=0.1時數(shù)值解結(jié)果如圖4和圖5： 圖4：H=0.1時孤子對在飽和介質(zhì)中的傳輸 圖5：H=0.1時孤子對傳輸?shù)慕孛鎴D公式（2）中H=0.2時數(shù)值解結(jié)果如圖6和圖7： 圖6：H=0.2時孤子對在飽和介質(zhì)中的傳輸 圖7：H=0.2時孤子對傳輸?shù)慕孛鎴D3、仿真結(jié)果分析 圖2明顯地反映了該孤子對在Kerr類型介質(zhì)中由于相互作用呈現(xiàn)周期性離合的特征，周期長度約為30，在整個傳輸過程中每個孤子在分離后的幅度沒有變化。圖3是圖2的截面圖，它反映了傳輸過程中孤子幅度的變

4、化。然后我們讓H=0.1，方程進而變?yōu)樵诮橘|(zhì)折射率非線性飽和條件下的非線性薛定諤方程。圖4是計算結(jié)果，可以看出在介質(zhì)折射率非線性飽和條件下孤子對的傳輸仍呈現(xiàn)周期性離合的特征，但周期長度明顯減小，約為17。圖5是圖4的截面圖，和圖3比較可知，兩孤子的匯合處不僅幅度增大而且范圍變寬。當(dāng)我們進一步增大介質(zhì)非線性飽和效應(yīng)時，令H=0.2,結(jié)果如圖6和圖7所示，與圖4和圖5比較，盡管周期長度和孤子幅度沒有太大的變化，但兩孤子相互作用的強度和范圍變得更大，且孤子波形的歧變也日趨嚴(yán)重。4、MATLAB程序 clearclc%h是空間間隔，nt是時間間隔，采用wspace函數(shù)把坐標(biāo)有時域變到頻域%A是H=0時

5、NLS方程的數(shù)值解%A1是H=0.1時NLS方程的數(shù)值解%A1是H=0.2時NLS方程的數(shù)值解h=0.01; z=0:h:80; T = 40; nt = 210; dt = T/nt; t = (1:nt)'-(nt+1)/2)*dt; w = wspace(T,nt); w=w' gn=-i*w.2;%以下程序計算介質(zhì)非線性飽和因子H=0時的能量分布A(1,:)=sech(t+3)+sech(t-3); H=0;for k=1:1:length(z)-1 A(k+1,:)=ifft(exp(gn*h/2).*fft(A(k,:).*exp(i*h*(abs(A(k,:).2

6、-H*abs(A(k,:).4);end%以下程序計算介質(zhì)非線性飽和因子H=0.1時的能量分布h=0.005; z1=0:h:40;A1(1,:)=sech(t+3)+sech(t-3); H=0.1;for k=1:1:length(z1)-1 A1(k+1,:)=ifft(exp(gn*h/2).*fft(A1(k,:).*exp(i*h*(abs(A1(k,:).2-H*abs(A1(k,:).4);end%以下程序計算介質(zhì)非線性飽和因子H=0.2時的能量分布A2(1,:)=sech(t+3)+sech(t-3); h=0.005; z2=0:h:40;H=0.2;for k=1:1:l

7、ength(z2)-1 A2(k+1,:)=ifft(exp(gn*h/2).*fft(A2(k,:).*exp(i*h*(abs(A2(k,:).2-H*abs(A2(k,:).4);endz=z'savefigureplot(t(1:1:length(t),A(1,:)xlabel('t')ylabel('A')title('光脈沖A=sech(t+3)+sech(t-3)')figurewaterfall(t(1:length(t),z(1:(length(z)-1)/50:length(z)-1),abs(A(1:(length(

8、z)-1)/50:length(z)-1,:)xlabel('t')ylabel('x')zlabel('|q|')title('H=0時孤子對在飽和介質(zhì)中的傳輸')view(0,45)figurestem(0:2:79,abs(A(1:(length(z)-1)/40:length(z)-1,526)xlabel('x')ylabel('|q|')title('H=0時孤子對在飽和介質(zhì)中傳輸?shù)慕孛鎴D')figurez=z1'waterfall(t(1:length(t),z

9、(1:(length(z)-1)/50:length(z)-1),abs(A1(1:(length(z)-1)/50:length(z)-1,:)xlabel('t')ylabel('x')zlabel('|q|')title('H=0.1時孤子對在飽和介質(zhì)中的傳輸')view(0,45)figurestem(0:1:39,abs(A1(1:(length(z)-1)/40:length(z)-1,515)xlabel('x')ylabel('|q|')title('H=0.1時孤子對在飽和

10、介質(zhì)中傳輸?shù)慕孛鎴D')figurez=z2'waterfall(t(1:length(t),z(1:(length(z)-1)/50:length(z)-1),abs(A2(1:(length(z)-1)/50:length(z)-1,:)xlabel('t')ylabel('x')zlabel('|q|')title('H=0.2時孤子對在飽和介質(zhì)中的傳輸')view(0,45)figurestem(0:1:39,abs(A(1:(length(z)-1)/40:length(z)-1,515)xlabel(&#

11、39;x')ylabel('|q|')title('H=0.2時孤子對在飽和介質(zhì)中傳輸?shù)慕孛鎴D')以下是文中用到的wspace函數(shù)的程序：function w = wspace(t,nt); % This function constructs a linearly-spaced vector of angular % frequencies that correspond to the points in an FFT spectrum. % The second half of the vector is aliased to negative %

12、frequencies. % % USAGE % % w = wspace(tv); % w = wspace(t,nt); % % INPUT % % tv - vector of linearly-spaced time values % t - scalar representing the periodicity of the time sequence % nt - Number of points in time sequence % (should only be provided if first argument is scalar) % % OUTPUT % % w - v

13、ector of angular frequencies % % EXAMPLE % % t = linspace(-10,10,2048)' % set up time vector % x = exp(-t.2); % construct time sequence % w = wspace(t); % construct w vector % Xhat = fft(x); % calculate spectrum % plot(w,abs(Xhat) % plot spectrum % % AUTHOR: Thomas E. Murphy () if (nar

14、gin<2) nt = length(t); dt = t(2) - t(1); t = t(nt) - t(1) + dt; end if (nargin = 2) dt = t/nt; end w = 2*pi*(0:nt-1)'/t; kv = find(w >= pi/dt); w(kv) = w(kv) - 2*pi/dt;參考文獻:1 方云團，王永順，沈廷根，李利方. 孤子對在飽和介質(zhì)中的傳輸J. 量子電子學(xué)報.2003，Vol.20 No.6.02 P.阿戈沃. 非線性光纖光學(xué)原理及應(yīng)用M. 北京: 電子工業(yè)出版社,2002.12.3 C.查布拉,P.卡納爾. 工程數(shù)值方法(第5版)M. 北京: 清華大學(xué)出版社,2007.12.4 張曉光等. 非線性薛定諤方程數(shù)值解法中傅立葉正逆變換選取的討論J. 計算物理. 2003,Vol.20,No.3:267-272. 5 呂理想,張曉萍. 不同形式非