孤立波的matlab-的数值计算和模拟

1前言孤立波简称孤波，是一种特殊的波，又是不难在自然界中出现的波，具有保持其波形和速度不变的特点。能发生强烈的相互作用，但相互作用后仍能保持其各自特点、形状、速度不变或只有一些位相改变的那些孤波称为孤子。孤波被称为自然界的相干结构，混沌运动所呈现的是非线性系统中奇妙的无序状态，相干结构则反映了非线性系统中的惊人有序性，因而孤立子理论的产生与发展是非线性偏微分方程研究中的一个重要组成部分。虽然随着人们研究的深入，在通信领域，采用光孤子技术大大提高了通信容量和通信距离，但作为非线性的两个重要的不同部分，混沌确为更多人所熟知。这可能是因为孤立子的研究和学习需要更高的数学和物理知识，对数学物理方法的要求也比较高。但本文将以matlab数值模拟为主要方法，简单形象得讲述KdV方程的性质和特点，并以此论证理论推导的结果。2KdV方程的数值模拟2.1孤波运动的模拟1895年荷兰数学家科特韦格(D.J.Korteweg)和德弗里斯(G.de.Vries)在浅水和小振幅假定下得到了浅水波单向运动方程。KdV方程经过变换，可得到各中不同的形式。为了方便，本文取以下形式：;06=++xxxttuuuu此方程的孤波解为*：)(()⎥⎥⎦⎤⎢⎢⎣⎡−−=022sec2,xctxchctxu；因为KdV方程的行波解都为同向运动的，则在进行单个孤子子运动的模拟时，我们把KdV方程的初始条件设为：；)((xhxu2sec20,=)利用差分法：孤立波的matlab的数值计算和模拟Madeby03级饶泽浪and吴辛烨;222;2;232112'''211''11'hfffffhffffhfffΔ−+−=Δ+−=Δ+=−−−−在matlab6.5下编程模拟，程序1如下：dt=0.0001;dx=0.1;c=dt/dx^3;d=dt/dx;x=-6:0.1:25;n=length(x);y(:,1)=2*(sech(x')).^2;%初始条件h=plot(x(1:end),y(1:end,1),'linewidth',2);axis([-6,25,0,3])set(h,'EraseMode','xor','MarkerSize',18);y(3:(n-2),2)=-c/2*(y(5:n,1)-2*y(4:(n-1),1)+…2*y(2:(n-3),1)-y(1:(n-4),1))+y(3:(n-2),1)...-6*y(3:(n-2),1).*d.*(y(4:(n-1),1)-y(2:(n-3),1));forj=1:100000;y(3:(n-2),3)=-c*(y(5:n,2)-2*y(4:(n-1),2)+…2*y(2:(n-3),2)-y(1:(n-4),2))+y(3:(n-2),1)...-6*y(3:(n-2),2).*d.*(y(4:(n-1),2)-y(2:(n-3),2));y(3:(n-2),1)=y(3:(n-2),2);y(3:(n-2),2)=y(3:(n-2),3);ifmod(j,40)==0;set(h,'XData',x,'YData',y(:,2));drawnow;endend其运行结果如图一；改变初始值，令)(()xhxu3sec60,2=；运行结果如图二；图一图二两图比较可发现，v1v2,即越大，移动速度越快，波高越大，但宽度越窄。这和行波解表达的物理意义是相符的：和是由初始条件决定的常数,为孤波的初始位置，c为孤波解的传播速度，波高为c/2，波宽反比于cc0x0xc。由计算模拟得到的结果我们不难发现这些结论的正确性。下面我们进一步讨论孤立波形状保持不变的物理内涵和物理意义。2.2波形的色散方程和色散关系色散方程的形式各样，为了对应KdV方程，本文取对应色散方程为：;0=+xxxtuu我们同样利用matlab进行模拟（程序附后），为了和KdV的方程结果有所比较，我们仍取初始条件为；)((xhxu2sec20,=))运行结果如图三。由图我们不难发现初始时刻出现的波形会随着时间的推移发生变化，发生弥散。其体现的物理内涵为：色散方程的通解*为:)([](wtkxiAtxu−∑=exp,;色散方程对应的色散关系*为：;03=+kw则;2kkwv−==我们不难发现对不同的k值，其传播速度v不等。则初始时刻由不同分波合成的波形，因分波的传播速度不同，波形发生变化，发生弥散，最后消失。2.3波动中的非线性汇聚我们取方程,;06=+ttuuu推导可见*。利用matlab进行模拟（程序在后），仍取初始条件为)(()xhxu2sec20,=；可得到图四：由图不难发现波速随着振幅的增加而增大，所以随着波的传播，波包的上端会越来越陡，最后坍塌，这就意味着某种汇聚效应。观察汇聚效应和色散效应的模拟图，不难发现，色散效应和汇聚效应的传播速度相反。这就是孤立波能稳定传播的原因，一个线性的波动在传播时由于存在色散效应，波动不能持续。只有当一个线性波动中存在非线性的汇聚效应时，且只有这两种作用达到某种平衡时，才能出现波形稳定的孤立波。2.4利用数值模拟验证KdV的初值问题KdV方程在不同的初始条件下，孤子的数目和运动情况都有所不同，一般情况下是采用反散色方法进行解决*，但本文将采用数值模拟的方式直观的得到一些结论。图三图四在程序1中，取)(()(xhxhVxu22sec2sec0,==)时，得到了图一，得到的是1个稳定的孤子。现取)(()()xhxhVxu22sec6sec0,==时，相应的得到图五可以发现，在此初始条件下，给出了2个孤子，而;326;212×=×=则我们可猜想当V，即，为正常数。在此初始条件下，将得到个孤立子。为了验证，我们取(1+=NN)xhNNxu2sec10,+∗=N)(()()N)(()()()xhxhNNxu22sec4*3sec10,=+∗=，得到图五，不难发现这猜想的正确性，限于篇幅，读者可以自行带入其它数值，即可发现结果确是如此。现在取V即，得到图六。不难发现未产生孤波，只有色散波。当,1−=)((xhxu2sec0,−=),0V读者取其它值仍会得到同样的结果，则当,V不能产生稳定的孤立子。0用同样的方法，读者代入不同的初始值，则会发现当V,0()()11+−NNVNN时，初始条件在KdV方程下将给出个孤立子，但也有色散波，限于篇幅，在此我们就不做演示，读者可自行演示。N由图形我们得到的这些结论，与由反散射方法理论推导的结果是一致的。2.5孤立子在KdV方程下的相互作用KdV方程只能实现单向运动，速度不同且互相追赶的多孤子碰撞。为了达到这个条件，我们依据孤波解的物理意义，取方程一的初始条件为：)(()())4(sec2)12(3sec60,22+++=xhxhxu；则可得到两孤立子的碰撞模拟，如图七；再取初始条件为：)(()++=)21(3sec60,2xhxu()())8(sec2)14(2sec422+++xhxh得到图八，其对应的就是三孤立子碰撞的情况。由图七，图八，我们不难发现：孤立子在碰撞的时候不满足一般线性波动的叠加原理。一般线性波动叠加遵循振幅叠加原理，但孤立子的碰撞，由图可发现振幅非但没有‘加’，反而‘减’了。碰撞过程就像大的孤子把小的孤子吞掉后，然后又把小的孤子吐了出来，并且各自都毫发无伤。这种现象很显然的是一种非线性的叠加，这也正是孤立波最重要的性质之一。至此我们把孤立子在一维KdV方程下的形成过程，及其在KdV方程下的一些最重要的性质，完全以matlab数值模拟的方式做了形象且正确的介绍。但是到此为止我们都是在一维的KdV方程下对孤立波进行模拟。下面我们将分别模拟在二维的KdV方程（即K-P方程）及Boussinesq方程下的孤波运动。图四（1,2,3）图五（1,2,3）图六（1,2,3）图七（1,2,3）3K-P方程的数值模拟作为二维的KdV方程，和KdV方程一样，K-P方程有着多种不同形式，但以两种形式为主。3.1KPI方程的数值计算及模拟为了和本文的KdV方程保持一致。我们取KPI方程为：()036=−++yyxxxxxtUUUUU；我们取它的初始条件*为：2220220220220]/1)(16)(4[/1)(16)(416)0,,(iiiikyykxxkyykxxyxu+−+−+−+−−=运行程序，得到相应的图九：由图我们不难发现，孤立子在KPI方程可以出现波包的形状，且碰撞时也是非线性叠加。但和KdV方程不同，它对横向的扰动不具有稳定性，具体可见*。3.2KPⅡ方程的模拟KPⅡ与KdV方程一样，它拥有孤波解：()[{]}0222234sec2),,(δλλ++−+=tkyxkhktyxU双孤立子的解为：),,(ln2),,(22tyxFxtyxU∂∂=；++=)1exp(1),,(ηtyxF)1221exp()2exp(A+++ηηη；()[]iiiiiitkyxkδλλη++−+=22342；()()()()22122122122144)12exp(λλλλ−−+−−−=kkkkA；利用孤波的解析解，我们进行模拟得到图十：利用双孤立子的解析解，我们进行模拟得到图十一：由图不难发现，KPⅡ不能产生波包，但是它有着与KdV很接近的性质，对横向扰动的稳定性*。图十图九图十一%线性色散图的差分解法clearclfdt=0.0001;dx=0.1;c=dt/dx^3;x=-100:0.1:100;n=length(x);y(:,1)=2*(sech(x')).^2;plot([-20,10],[0,0],'r');text(0,0,'水平面');holdonh=plot(x,y(:,1),'linewidth',2);axis([-20,6,-1,3])set(h,'EraseMode','xor','MarkerSize',18);y(3:(n-2),2)=-c/2*(y(5:n,1)-2*y(4:(n-1),1)+2*y(2:(n-3),1)-y(1:(n-4),1))+y(3:(n-2),1);forj=3:30000;y(3:(n-2),3)=-c*(y(5:n,2)-2*y(4:(n-1),2)+2*y(2:(n-3),2)-y(1:(n-4),2))+y(3:(n-2),1);y(3:(n-2),1)=y(3:(n-2),2);y(3:(n-2),2)=y(3:(n-2),3);ifmod(j,40)==0;set(h,'XData',x,'YData',y(:,2));drawnow;pauseendend%非线性汇聚的差分解法clearclflbd=-6;dx=0.1;dt=0.0001;x_begin=-4;x_end=4;x=[x_begin:dx:x_end];xlength=length(x);u(:,1)=2*(sech(x')).^2;u(2:(xlength-1),2)=lbd.*(dt/(2*dx)).*u(2:(xlength-1),1).*(u(3:xlength,1)-u(1:(xlength-2),1))+u(2:(xlength-1),1);h1=plot(x,u(:,1),'r','linewidth',2);axis([x_beginx_end-0.22.2])set(h1,'ydata',u(:,2));forj=3:850u(2:(xlength-1),3)=lbd.*(dt/dx).*u(2:(xlength-1),2).*(u(3:xlength,2)-u(1:(xlength-2),2))+u(2:(xlength-1),1);u(:,1)=u(:,2);u(:,2)=u(:,3);ifmod(j,400)==0;set(h1,'xdata',x,'ydata',u(:,3));drawnowendend%kp1的单孤波（行波法）clearc0=6;beta=6;dx=0.2;dy=0.2;dt=0.002;x_begin=-22;x_end=15;y_begin=-30;y_end=20;x=[x_begin: