1 / 34

Zhejiang Normal University

Zhejiang Normal University. Split-Step Fourier Transform Method. 李画眉. 2009 年 11 月. Outline. 一、 Algorithms 二、 Stationary solution for 1+1D NLSE 三、 Propagation for 1+1D NLSE 四、 Stability Analysis 五、 Conclusions. Algorithms ： Image-time method.

adam-monroe
Download Presentation

Zhejiang Normal University

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Zhejiang Normal University Split-Step Fourier Transform Method 李画眉 2009年11月

  2. Outline 一、 Algorithms 二、 Stationary solution for 1+1D NLSE 三、Propagation for 1+1D NLSE 四、 Stability Analysis 五、Conclusions

  3. Algorithms：Image-time method By making use of a substitution t-it, then one can employ other appropriate algorithms (FFT, difference scheme, etc.) to obtain stationary solution. Advantage: independence of initial guess value Disadvantage: only ground state can be obtained

  4. linear nonlinear linear Stationary solution for 1+1D NLSE • FFT Linear part Nonlinear part

  5. Program clear all p=1; omega=2; %%------------------- n=2048; hx=0.06; x=(-n/2:n/2-1)*hx; hw=2*pi/(n*hx); w=fftshift((-n/2:n/2-1)*hw); %%------------------- % q=exp(-1*(x).^2/2); q=sech(x); intensity=4.6; u1(:,1)=(abs(q).^2)'; %------------------- V=cos(omega*x);

  6. Program %-------------------- L=50; nm=5000; h=L/nm; %------------------- for j=1:nm j D=exp(((i*w).^2/2)*h/2); qstep1=ifft(D.*fft(q)); N=exp((p*V+(abs(qstep1)).^2)*h); qstep2=N.*qstep1; q=ifft(D.*fft(qstep2)); q=sqrt(intensity)*q/sqrt(sum(abs(q).^2)*hx); u=abs(q); r=floor(2+(j-1)/50); u1(:,r)=u'; end

  7. Program kin=-sum((q(3:end)-q(1:end-2)).^2)/4/hx; p_i=sum(2*q.^2.*(abs(q).^2+V))*hx; b=(kin+p_i)/2/intensity z=0:h*50:50; figure(1) mesh(x,z,u1'); view(0,90) figure(2) plot(x,u1(:,end),'r',x,V,'b')

  8. Energy_b Energy_b

  9. Program clear all p=1; omega=2; %%------------------- n=2048; hx=0.06; x=(-n/2:n/2-1)*hx; hw=2*pi/(n*hx); w=fftshift((-n/2:n/2-1)*hw); %%------------------- loop=0;

  10. Program for intensity=1:0.1:5 % q=exp(-1*(x).^2/2); q=sech(x); u1(:,1)=(abs(q).^2)'; %------------------- V=cos(omega*x); %-------------------- L=50; nm=1000; h=L/nm; loop=loop+1 %-------------------

  11. Program for j=1:nm j; D=exp(((i*w).^2/2)*h/2); qstep1=ifft(D.*fft(q)); N=exp((p*V+(abs(qstep1)).^2)*h); qstep2=N.*qstep1; q=ifft(D.*fft(qstep2)); q=sqrt(intensity)*q/sqrt(sum(abs(q).^2)*hx); end kin=-sum((q(3:end)-q(1:end-2)).^2)/4/hx; p_i=sum(2*q.^2.*(q.^2+p*V))*hx; b(loop)=(kin+p_i)/2/intensity end intensity=1:0.1:5 plot(b,intensity)

  12. Numerical results Iterative process Energy_b Stationary solution

  13. Propagation for 1+1D NLSE NLSE

  14. Program clear all p=1; omega=2; %%------------------- n=2048; hx=0.06; x=(-n/2:n/2-1)*hx; hw=2*pi/(n*hx); w=fftshift((-n/2:n/2-1)*hw); %%------------------- % q=exp(-1*(x).^2/2); q=sech(x); intensity=1; %------------------- V=cos(omega*x);

  15. Program %%%%%%%%%%%%%%%% L=50; nm=5000; h=L/nm; u1(:,1)=(abs(q).^2)'; %------------------- for j=1:nm j D=exp((i*(i*w).^2/2)*h/2); qstep1=ifft(D.*fft(q)); N=exp((i*p*V+i*(abs(qstep1)).^2)*h); qstep2=N.*qstep1; q=ifft(D.*fft(qstep2)); u=abs(q).^2; r=floor(2+(j-1)/50); u1(:,r)=u'; end

  16. Program z=0:50*h:L; mesh(x,z,u1'); view(0,90)

  17. Numerical results Propagation

  18. Stability Analysis for 1+1D NLSE Assuming the stationary solution is of the form We consider the stability of the stationary solution w(x) by employing the following algorithms: Eigenvalue method Split-step Fourier method

  19. Eigenvalue method In eigenvalue method, there are two assumptions for the perturbed stationary solution First assumption: Physica D 237, 3252 (2008) where the perturbation components u, v could grow with a complex rate λduring propagation. Substitute the perturbed solution into equation, we obtained the linear eigen-equations The solution is stable if the imaginary parts of λ equal zero.

  20. We begin by discretizing the domain x∈[-a,a] by placing a grid over the domain. We will use the grid with grid spacing h=2a/N in axis x. We will attempt to approximate the solution at the points on this lattice, and define uk and vk to be functions defined at the point xk≡-a+(k-1)h or the lattice point k, where k=1, 2, ⋯, N+1. Thus we can obtain the difference scheme as follows with Espically and due to

  21. namely where The soliton is stable if the imaginary parts λ of equal zero. stability_1D_eigenvalue_1.m

  22. Second assumption: PRL93, 153903 (2004) with the complex rate δ: After substituting this perturbed solution into equation and linearizing, we obtain a linear eigenvalue problem as follows The soliton is stable if the real parts of δequal zero.

  23. Discretizing Namely where stability_1D_eigenvalue_2.m

  24. Discretizing Namely where stability_1D_eigenvalue_2.m

  25. Split-step Fourier method JOSAB21, 973 (2004) By denoting where w(x) is the fundamental soliton and U(x,z)<<1 is the infinitesimal perturbation, the linearized equation for U(x,z)is Starting from a white-noise initial condition, we simulate above linear equation for a long distance (hundreds of z units).

  26. FFT RK method FFT Note FFT RK method

  27. Program clear all tic p=1;omega=2; %%------------------- n=1024;hx=0.08; x=(-n/2:n/2-1)*hx; hw=2*pi/(n*hx); w=fftshift((-n/2:n/2-1)*hw); %%------------------- q=exp(-1*(x).^2/2); q=sech(x); intensity=4.6; u1(:,1)=(abs(q).^2)'; %------------------- V=cos(omega*x);

  28. Program %-------------------- L=50; nm=5000; h=L/nm; %------------------- for j=1:nm j D=exp(((i*w).^2/2)*h/2); qstep1=ifft(D.*fft(q)); N=exp((p*V+(abs(qstep1)).^2)*h); qstep2=N.*qstep1; q=ifft(D.*fft(qstep2)); q=sqrt(intensity)*q/sqrt(sum(abs(q).^2)*hx); u=abs(q); r=floor(2+(j-1)/50); u1(:,r)=u'; end

  29. Program kin=-sum((q(3:end)-q(1:end-2)).^2)/4/hx; p_i=sum(2*q.^2.*(abs(q).^2+V))*hx; b=(kin+p_i)/2/intensity %%-------------------------- alpha=1/(2*hx^2); w=u1(:,end)'; beta=2*w.^2+p*cos(omega*x)-b-1/hx^2; kai=w.^2; growth_rate1=stability_1D_eigenvalue_1(alpha,beta,kai,w,n) %%------------------- % alpha=1/(2*h^2); % w=phi; p=2*alpha-p*cos(omega*x)+b; % kai=w.^2; growth_rate2=stability_1D_eigenvalue_2(alpha,p,kai,w,n) toc

  30. Subprogram1 function growth_rate=stability_1D_eigenvalue_1(alpha,beta,kai,phi,N) A=[beta(2) kai(2);-kai(2) -beta(2)]; B=[alpha 0;0 -alpha]; for kk=3:N A=[A [zeros(2*(kk-3),2);B];zeros(2,2*(kk-3)) B [beta(kk) kai(kk);-kai(kk) -beta(kk)]]; end growth_rate=max(abs(imag(eig(A))));

  31. Subprogram2 function growth_rate=stability_1D_eigenvalue_2(alpha,p,kai,phi,N) for mm=1:N-2 AA2(mm,mm)=-p(mm+1)+3*kai(mm+1); AA2(mm,mm+1)=alpha; AA2(mm+1,mm)=alpha; BB2(mm,mm)=p(mm+1)-kai(mm+1); BB2(mm,mm+1)=-alpha; BB2(mm+1,mm)=-alpha; end AA2(N-1,N-1)=-p(N)+3*kai(N); BB2(N-1,N-1)=p(N)-kai(N); CC2=[zeros(N-1,N-1) AA2;BB2 zeros(N-1,N-1)]; delta2=eig(CC2); NN2=length(delta2); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% growth_rate=max(abs(real(delta2)));

  32. Numerical results stability_eigen_newton_1D_LU_growth.m

  33. Conclusions • As far as ground state is concerned, image-time method is a most effective algorithm to deal with it due to insignificance of initial guess value. It is worth noting that energy conservation law must be met when image-time method is employed. • Split-step Fourior algorithm is a feasible method for stability analysis.

  34. 谢谢大家！

More Related