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

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