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

## # 牛顿冷却定理

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

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

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

## # 一维波动方程

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

$u(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

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

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

## # 二维波动方程

$\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)$

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

$u(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

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

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
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

[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)