# 波动方程的数值解(Mathematica)

今日我尝试使用神经网络去解微分方程。为了去验证神经网络解法的正确性,我想要将神经网络的结果与数值 解进行对比。但是,对解微分方程的数值解只会用欧拉法的我来说,解稍微复杂一些的微分方程就特别困难。 网络上也有很优秀的使用Python求解微分方程数值解的方法,但是我很懒,懒得自己写程序了。 由于我的目的并不是想要学习微分方程的数值解法,而是仅仅用数值解与神经网络的数值解进行对照,因此我 选择使用Mathematica来作为我需要的工具。Mathematica给我的感觉是面向问题的处理工具,可以让我 专注于我要解决的问题,帮助我快速实现,验证我的想法。

# 牛顿冷却定理

为了熟悉Mathematica的数值解工具NDSolve,我们先试一个经典的微分方程: 牛顿冷却定律 (Newton's law of cooling)。这个冷却 定理说的是:

一个物体和其周围处于一个不同的温度下的话,最终这个物体会和其周围达成一个相同的温度。一个比较 热的物体将会冷却,因为它使其周围变温暖。一个比较冷的物体会因为其周围的高温而温度上升。

微分方程的描述是下面这个样子:

dTdt=k(TC)\frac{dT}{dt} = -k(T-C)

其中T代表物体的温度,t是时间,k是一个和物体表面积和传热系数相关的一个常数,c代表的是 环境的温度。这个微分方程的意思可以这样通俗地解释为: 一个物体温度下降的速率dTdt-\frac{dT}{dt} 和这个物体和环境的温度差TCT-C的大小成正比。至于有些人可能马上就会说的,你这不对啊。有些时候把 开水和常温水放进冰箱,会出现开水先结冰这样的现象。对的,这就是传说中的姆潘巴现象 (Mpemba effect)。 可见这个微分方程是有一定的局限性的。但这不是我们现在要讨论的问题。

现在我们假设环境的温度是c=25°C, 物体的初始温度是T(0)=80°C, 常数k=2现在我们想要推测出 这个物体的温度随时间的变化关系。也就是找到函数不同时间t物体的温度T(t)。 因此微分方程的解是一个函数。

现在开始Mathematica代码:

k := 2; 
c := 20;
s = NDSolve[{D[T[t], t] == -k*(T[t] - c), T[0] == 80}, T, {t, 0, 10}];
Plot[Evaluate[{T[t] /. s}], {t, 0, 10}, PlotRange -> {{0, 10}, {10, 80}}, AxesLabel -> {t, T}]
1
2
3
4

这里的NDSolve函数是得到微分方程数值解的关键。Mathematica中的函数表示不同与Python用括号表示函数。 Mathematica中的函数用中括号[]表示。主要是要告知待解的微分方程D[T[t], t] == -k*(T[t] - c), 初始条件 T[0] == 80。数值解T的范围告知Mathetatica, {t, 0, 10}。 最后将计算的结果放在s中。接下来一行是作图,其中Evaluate[{T[t] /. s}]代表获得s中的数值。

图1:物体起始时刻的温度为80°C条件下微分方程的解

改变不同的物体初始温度,我们可以得到下面的结果:

k := 2; 
c := 25;
s0 = NDSolve[{D[T[t], t] == -k*(T[t] - c), T[0] == 80}, T, {t, 0, 10}];
s1 = NDSolve[{D[T[t], t] == -k*(T[t] - c), T[0] == 40}, T, {t, 0, 10}];
s2 = NDSolve[{D[T[t], t] == -k*(T[t] - c), T[0] == 20}, T, {t, 0, 10}];
s3 = NDSolve[{D[T[t], t] == -k*(T[t] - c), T[0] == 10}, T, {t, 0, 10}];
img = Plot[{Evaluate[{T[t] /. s0}], Evaluate[{T[t] /. s1}], 
   Evaluate[{T[t] /. s2}], Evaluate[{T[t] /. s3}]}, {t, 0, 10}, 
  PlotRange -> {{0, 10}, {10, 80}}, AxesLabel -> {t, T}]
1
2
3
4
5
6
7
8
9
Export["~/Desktop/img.svg", img]
1
图2:不同初始条件下对应的微分方程的解

可以看到,当物体的温度T高于环境温度c=25时,T随时间减小,最后稳定在环境温度;反之,温度升高后趋于稳定。

# 一维波动方程

振动的吉他弦可以近似地用一维波动方程来描述:

2ut2=c22ux2\frac{\partial^{2} u}{\partial t^{2}}=c^{2} \frac{\partial^{2} u}{\partial x^{2}}

这里的u表示的是吉他弦偏离平衡位置的距离,t是时间,x指琴弦上的点到上弦枕的距离。c代表弦上的振动波 传播的速度,c=T/μc = \sqrt{T / \mu}, 它由弦的张力T,和单位长度上弦的质量μ\mu两个因数共同决定。

要解这个微分方程,还需要告知初始时刻的位置和速度,以及边界条件。首先初始位置条件我们暂时这样表示:

u(0,x)=0.5e100(x1)2u(0, x) = 0.5*e^{-100*(x - 1)^2}

img = Plot[0.5*E^(-100*((x - 1)^2)), {x, 0, 2}, PlotRange -> {{0, 2}, {-0.5, 0.5}}, AxesLabel -> {x, u}];
Export["~/Desktop/init0.svg", img]
1
2
图3: 初始位置

初始时刻弦上所有点的速度为0,因此有utt=0=0\frac{\partial u}{\partial t}|_{t=0} = 0。边界条件,即弦两端偏离平衡 位置的位移为0,表示为

u(t,0)=u(t,2)=0u(t,0) = u(t,2) = 0

我们这个例子中认为弦长为2;接下来,我们将主方程,初始条件,边界条件带入,使用Mathematica求解:

s = NDSolve[{D[u[t, x], t, t] == 0.02*(D[u[t, x], x, x]), 
    u[t, 0] == 0, u[t, 2] == 0, u[0, x] == 0.5*E^(-100*((x - 1)^2)), 
    Derivative[1, 0][u][0, x] == 0}, u, {t, 0, 20}, {x, 0, 2}];
gif = Table[
   Plot[Evaluate[u[t, x] /. s], {x, 0, 2}, AxesLabel -> {x, u}, 
    PlotRange -> {{0, 2}, {-0.5, 0.5}}, 
    PlotLabel -> "t=" <> ToString[NumberForm[t, {2, 2}]]], {t, 0, 20, 
    0.5}];
Export["~/Desktop/1d_new.gif", gif]
1
2
3
4
5
6
7
8
9

于是我们得到微分方程的解

图4: 一维波动方程的解

从这个解中我们可以清楚地看到波的传播,反射,和干涉。这里需要提的是,虽然我们一直拿吉他弦的振动来 形象地展示波动方程,但是这里的解其实并不代表实际中吉他弦的振动方式。不同的点在于我们这里给的初始 位置条件不同于我们拨弦时的初始位置条件。上面例子中用到的初始位置条件类似一个“钟形”,而实际 当中拨吉他弦类似于一个“三角形”,拨弦手指的位置位于三角形的顶点位置。那么这样的初始位置条件必然就 应该是一个分段函数。

# 二维波动方程

接下来,我们从一维波动方程扩展到二维。二维波动方程可以用来描述扔石子到水中,或者是鼓面的振动。

2ut2=c2(2ux2+2uy2)\frac{\partial^{2} u}{\partial t^{2}}=c^{2}\left(\frac{\partial^{2} u}{\partial x^{2}}+\frac{\partial^{2} u}{\partial y^{2}}\right)

二维波动方程中,点(x,y)沿竖直方向z的位移表示为u(t,x,y)u(t,x,y), 我们限制x[0,2],y[0,2]x\in[0,2], y\in[0,2], 因此, 边界条件:

u(t,0,y)=u(t,2,y)=u(t,x,0)=u(t,x,2)=0u(t,0,y) = u(t,2,y) = u(t,x,0) = u(t,x,2) = 0

初始时,我们给定一个初始位置,模拟戳了一下水平面,仿照上面一维波动方程,给定下面的初始位置条件:

u(0,x,y)=0.5e100(x1)2100(y1)2u(0, x, y) = 0.5*e^{-100*(x - 1)^2 -100*(y - 1)^2}

img = Plot3D[0.5*E^(-100*((x - 1)^2)-100*((y - 1)^2)), {x, 0, 2}, {y, 0, 2}, PlotRange -> {{0, 2},{0,2}, {-0.5, 0.5}}, AxesLabel -> {x, y, u}];
Export["~/Desktop/init2.svg", img]
1
2
图5: 2维波函数的初始位置

此外,我们认为初始速度处处为0,于是得到速度初始条件:

utt=0=0\frac{\partial u}{\partial t}|_{t=0} = 0

带入主方程和各个条件到NDSolve得到数值解:

s = NDSolve[{D[u[t, x, y], t, t] == 
    0.02*(D[u[t, x, y], x, x] + D[u[t, x, y], y, y]), u[t, 0, y] == 0,
    u[t, 2, y] == 0, u[t, x, 0] == 0, u[t, x, 2] == 0, 
   u[0, x, y] == 0.5*E^(-100*((x - 1)^2 + (y - 1)^2)), 
   Derivative[1, 0, 0][u][0, x, y] == 0}, 
  u, {t, 0, 12}, {x, 0, 2}, {y, 0, 2}]
1
2
3
4
5
6
gif = Table[
  Plot3D[Evaluate[u[t, x, y] /. s], {x, 0, 2}, {y, 0, 2}, 
   AxesLabel -> {x, y, u}, PlotRange -> {{0, 2}, {0, 2}, {-0.5, 0.5}},
    PlotLabel -> "t=" <> ToString[NumberForm[t, {2, 2}]]], {t, 0, 12, 
   0.5}];
Export["~/Desktop/2dWave.gif", gif]
1
2
3
4
5
6
图4: 二维波动方程的一个特解
s = NDSolve[{D[u[t, x, y], t, t] == 
     0.02*(D[u[t, x, y], x, x] + D[u[t, x, y], y, y]), 
    u[t, 0, y] == 0, u[t, 2, y] == 0, u[t, x, 0] == 0, 
    u[t, x, 2] == 0, 
    u[0, x, y] == 
     0.5*E^(-100*((x - 0.5)^2 + (y - 0.5)^2)) + 
      0.3*E^(-100*((x - 1.5)^2 + (y - 1.5)^2)), 
    Derivative[1, 0, 0][u][0, x, y] == 0}, 
   u, {t, 0, 12}, {x, 0, 2}, {y, 0, 2}];
gif = Table[
   Plot3D[Evaluate[u[t, x, y] /. s], {x, 0, 2}, {y, 0, 2}, 
    AxesLabel -> {x, y, u}, 
    PlotRange -> {{0, 2}, {0, 2}, {-0.5, 0.5}}, 
    PlotLabel -> "t=" <> ToString[NumberForm[t, {2, 2}]]], {t, 0, 12, 
    0.5}];
Export["~/Desktop/Cool2dWave.gif", gif]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
图5: 修改初始条件后二维波动方程的另一个特解

最后, Mathematic Notebook 文件可以在这里 (opens new window)下载。

参考:

[1] 维基百科: 冷却定率 (opens new window)
[2] 维基百科: 姆潘巴现象 (opens new window)
[3] Giordano, N. (1998). Physics of vibrating strings. Computers in Physics, 12(2), 138–145. https://doi.org/10.1063/1.168621
[4] 知乎: Python数值计算:二维波动方程有限差分解 (opens new window)

上次更新: 11/24/2021, 10:39:29 PM