“孤子”及计算机数值方法

1““““孤子孤子孤子孤子””””及计算机数值方法初探及计算机数值方法初探及计算机数值方法初探及计算机数值方法初探刘景博2008011067清华大学无810班摘要“孤子”是非线性科学里一朵美丽的奇葩，被数学家和物理学家所津津乐道；而光孤子通信被认为是未来最有前景的通信技术之一。但由于孤子方程的非线性性，单纯用解析方法研究孤子往往很难。本文主要用数值方法展示了孤子的一些基本性质：首先计算了浅水波KdV方程的孤子解，然后通过数值实验验证了“孤子”的波形不变性、碰撞分离的稳定性；最后数值模拟了1阶、2阶和3阶孤立子光脉冲波形在光纤中的演化过程，并根据图形对结果进行了归纳总结。关键词：孤子、非线性科学、KdV方程、光孤子通信、NLS方程一、“孤子”简介1.“孤子”的发现1834年，苏格兰海军工程师罗素(J.ScottRussell)在沿着狭窄的运河骑马时发现了“孤波”现象。他观察到两匹骏马拉着一条船沿运河迅速前进，当船突然停止时，随船一起运动的船头处的水堆并没有停止下来，而是激烈地在船头翻动起来，随即突然离开船头，并以巨大的速度向前推进。一个轮廓清晰、光滑的水堆，犹如一个大鼓包，沿着运河一直向前推进，并且在行进过程中其形状与速度没有明显变化。罗素纵马追逐了两英里，它才渐渐消失。罗素当时将它称为“孤波”，现在一般称为“孤子”。罗素在水槽中反复地做了实验，但未能对“孤子”作出理论解，也没有说服自己的同事们。这个事件于是引起了断断续续的争论，直到1895年，两位荷兰科学家科特维格(Kortweg)与德弗雷斯(deVries)导出了浅水波的微分方程（KdV方程），并得出了一个类似孤子的解析解，对孤子的争论才渐渐平息[1]。2.为什么研究“孤子”？在KdV方程被导出后的几十年里，“孤子”似乎被人们遗忘，但它不会永远沉默。1955年，为了从数值实验上验证统计力学中的能量均分定理，著名物理学家费米(E.Fermi)、帕斯塔(J.Pasta)和乌莱姆数值计算了用非线性弹簧联结的64个质点组成的弦的振动。计算结果令人意外，长时间以后能量几乎全部又回到了少数质点上。由于非线性振子系统的能量不均分问题可以同KdV方程联系起来。2他们于是发展了一套数值和解析相结合研究非线性方程方法，即从数值结果和图形显示中获得定性启示，再尝试用解析方法给与证明，然后再用数值分析检验解析的推论，如此循环，步步深入。乌勒姆将这种方法称为“计算协同学”。在之后的短短几年里，人们获得了大量关于KdV方程和一大批非线性偏微分方程的解析结果。研究表明，包括KdV方程的许多偏微分方程都有孤立波解。可见孤立波是自然界里一种既特殊而又不难见到的波动现象。3.“孤立子”的性质和它在光通信中的应用1）碰撞稳定性科学家仍采用数值模拟的方法，发现两个孤立波可以相遇、碰撞，而且分离之后波包的形状不会发生明显的变化。人们把具有稳定碰撞特性的孤立波称为“孤立子”或称“孤子”(soliton)。孤立子光脉冲的碰撞稳定性为设计波分复用提供了方便。2）孤子脉冲的不变性由于孤立子光脉冲在光纤中传播时具有稳定不变的能量与波形，人们想到了将孤子应用于光纤通信技术。我们注意到，目前的光纤通信技术采用低强度光脉冲的线性通信方式。低强度光脉冲在光纤中传播不可避免地产生色散，从而造成光脉冲的加宽与变形，这大大地影响到光信息传送的质量与距离。为了长距离、高质量的传送信息，必须在传送路程上设置许多造价很昂贵的的中继站。采用孤立子进行通信为解决这一问题提供了新的思路。光学孤立子虽然在传播时能够保持稳定的能量和波形，但是仍有某些因素（例如光纤内部的微小瑕疵）可能造成孤立子的能量损失。不过人们想到了在光的传播过程中给孤立子补充能量的办法，而不用设置中继站。拉曼泵浦技术可以实现对孤立子能量损失的补充。该技术的实质是：当两列不同频率的光波在光纤中共同传输时，如果它们的能量足够大，高频的光波会将其部分能量转移给低频的光波，其基本工作原理如图1所示[2]。在传播过程中，光学孤立波将通过与泵浦波（980nm）的相互作用而获得能量。1988年，采用周期性增益补偿，进行了4000km的长距离试验，现在已有成功地进行了近万km长距离试验的报道。由于孤立子的种种奇特的性质，在未来长距离、高速、大容量通信中，光孤立子通信展现出了诱人的美妙前景。3二、KdV方程的孤子解的数值分析1.KdV方程的行波解由流体力学得到的KdV方程是描述浅水波的微分方程，形式如下：033=∂∂+∂∂+∂∂xuxuutu(1)上式中，u代表水波相对于静止水面的高度，即波幅，x和t分别为位置和时间坐标。KdV方程是一个非线性的偏微分方程，求解是很困难的。我们先考虑它的行波解。假设一个人以与孤子相同的速度追逐孤子，那么人会看到波面静止不变。也就是说，在惯性系S1(x1,t)中，u=u(x1,t)满足tu∂∂=0，其中s1与静止系0s有变换关系x1=x-vt。这样可以把原来的偏微分方程转化为常微分方程。令x1=x-vt，t=t。则有txxtxuttxuttxu∂∂∂∂+∂∂=∂∂1111),(),(),(=111),(),(xtxuvttxu∂∂−∂∂=∂∂xu1xu∂∂，31333xuxu∂∂=∂∂将他们带入原方程（1），并由tu∂∂=0可得(u-v)3131xuxu∂∂+∂∂=0上述方程可以写成常微分方程。两端同乘dx1，得0''=+−duvduudu信号光（1550nm）波分复合器泵浦光980nm光纤1500nm980nm1550nm980nm隔离器1550nm图1拉曼泵浦4积分得：Cuvuu=+−''212。显然当∞→1x时，0→u，0→′′u。于是有C=0，方程化为0''212=+−uvuu(2)给等式(2)两端同乘1'dxu，得0'''2112=+−dxuuvududuu两端积分得：cuvuu=+−223'212161。令1x∞→,由0→u，0'→u得到c=0.即方程化为0'212161223=+−uvuu由此解出3231'uvuu−±=，即cdxuvudu+±=−∫1323/于是我们得到行波解u(x1)=3vsech2(12xv)(3)下面介绍在Matlab使用龙格——库塔公式求解常微分方程(2)的方法。对于形如(2)的二阶常微分方程，先将它转化为一阶常微分方程组：+−==1211221121vuudxduudxdu(4)初值条件为：u(0)=A,0)0('=u。这样画出的图形应该是孤子的一半（即x0的部分）。在Matlab中用使用龙格——库塔公式计算常微分方程的函数调用格式为[3]：[x,y]=ode45(‘Fun’,Tspan,y0,options)，其中括号内的前三个参数分别为导数的表达式文件、自变量取值范围、函数初值。由前面的论述知，单个孤子解的表达式为：u(x1)=3vsech2(12xv)我们将这个理论结果和数值计算所得的微分方程（2）的解画在一张图上，做一5比较。取v=1，运行后得到图2。可见当x16时两条曲线几乎重合。改变v的值，两条曲线始终很接近。这说明数值解法在自变量变化范围不大的情况下有很高的精度。另外如果我们尝试改变u的初值，就会发现方程(2)其他的解曲线不一定都是“孤子”形状（图略），这表明孤子解只是方程(2)的一类特殊的解。024681000.511.522.53xuSolutionofODEforv=1数数数数解数图2方程(2)的数值解和解析解2.“孤子”波形不变和碰撞稳定性的验证我们考虑u(x,t)的数值解，这样便可以观察波运动的过程，看它是否真能保持波形不变，并且有碰撞分离的稳定性。为了用差分法描述微分方程(1)，我们先定义一个矩阵u(i,j)来表示u的值。假设x(i)和t(j)位置与时间的离散值，u(i,j)表示时间为t(j),在位置x(i)处的u值。如果给定初值，再让i,j跑遍某一区间内的值，求出对应的u(i,j)，我们就知道了孤子这一区间内的运动情况。假定),0(),,0(mmttxx∈∈，∆t和∆x分别为时间与位置的步长。在))(),((jtix处，一阶导数用中心差分代替，即∆−−+≈∂∂∆−−+≈∂∂)2/()),1(),1(()2/())1,()1,((xjiujiuxutjiujiutu(5)3阶导数33xu∂∂用多项式插值法来求。固定时间t，u是x的函数)(xu。过其上五点6(x1,u1),(x2,u2),(x3,u3),(x4,u4),(x5,u5)的差值多项式为∑∏≤≤≠−−51)()(kkikikkxxxxu，我们将该多项式在3x处的3阶导数作为)(xu在x3处的3阶导数。由此推出在x(i)，t(j)处有333)/()2/),2(),1(),1(2/),2((xjiujiujiujiuxu∆+++−−+−−≈∂∂(6)由(5)与(6)式得到方程（1）的如下离散显格式：)2/()),1(),1()(,((2)1,()1,(xjiujiujiutjiujiu∆−−+∆−−=+))/()2/),2(),1(),1(2/),2((3xjiujiujiujiu∆+++−−+−−+NjMi,,2,1,0,,,2,1,0LL==(7)基于离散格式（7）的方程（1）求解步骤如下：Step1给定定解条件。初始条件为u(x1，0)=3vsech2(12xv)。对于x=0,和x=mx的边值条件，由周期性延拓(u(x,t)=u(x+xm,t))而确定。Step2选取时间步长t∆与距离步长x∆。Step3用离散格式递推计算),(jiu。图3给出了利用上述算法得到的方程（1）的数值解，其中2.0=∆t，1=∆x，3.0=v。可以看到，随着时间t的增加，波形几乎没有任何改变，而且波包在tm=200的时间内大约运动了60个长度单位，这与程序中设定的v=0.3吻合得很好。这就“验证”了孤子的稳定性。下面我们看一看两个速度不同的水面孤子相碰的情况。由于孤子的局域性，当两个孤子相距很远时可看作两个独立的孤子，因此只需将初值改为两个相距较远、速度不同的孤子波形的叠加就可以了（当两个孤子相距很近时不能这样叠加）。我在上述算法中令v1=0.1,v2=0.5,即得图4，给出了两孤子的相遇过程。图5为图4中t=82，102，122，142，162的横截图。由图5可见，速度快的那个波更高、更窄。值得注意的是，在两个波相遇处水波的高度反而低于最高的那个波的高度！事实上这是由KdV方程的非线性所致。一般来说，在非线性系统中两个波相遇时不能简单的叠加。另外由图4还可以看到，两个孤立波相遇后其轨道相对于原来略有偏移。这说明碰撞导致了相移。我们看到了计算机数值计算的威力：孤子的一些性质已经在图上直观、清晰地表示出来。更重要的一点是，非线性方程是个性很强的问题。如果我把方程(1)的形式改变一点点，可能求解的方法和解的表达式变化很大，或根本不能求得解析解。数值计算似乎为我们提供了一种较为通用而且有效的方法，尽管在有些情况下它不能从根本上证明一些问题(例如为什孤立子碰撞后保持形状不变)。7图4两个孤子的相遇过程图3单个孤子的运动情况8图5两个孤子的相遇过程三、光孤子的性质的数值分析在这个例子中，我们会进一步体会到计算机在现代科学研究中不可替代的作用。我们知道光脉冲可以作为信息的载体，因此如果光纤中能存在光孤子，那么它就可以用来长距离地传输信息。人们发现，某种条件下光纤中的非线性效应会补偿能使脉冲展宽的群速色散效应，因而产生光孤子。光孤子遵循非线性薛定谔方程(NLSE)[6]i021222=+∂∂+∂∂uuuuτξ(8)为了方便起见，以下用x、t代替符号ξ、τ。人们已经通过理论推导得到了该方程的一阶孤立子解：2/)(sech),(ixettxu=,并且知道如果在x=0处输入形如u(t)=Asech(t)（A=N为整数）的脉冲就能得到光孤子。当N=1时，对应的一阶光孤子可以保持稳定的波形；当N1时孤子会以x=2π为周期发生波形的改变。不过即使是N=2对应二阶光孤子，它的解析表达式也是相当复杂的。这时数值计9算有派上了用场。下面就从(8