一维 Burgers 方程如下

\[ \left\{ \begin{aligned} u_t + \left(\frac{u^2}2\right)_x &= 0\\ u(x,0) &= u_0(x) \end{aligned} \right. \]

其中定义域 \(x \in \mathbb{R}\)，通常假定初值有周期性，精确解此时也具有周期性。

特征线方程

从\((\xi,0)\)发出的曲线记作 \(l_\xi\)，曲线参数为 \(s\)，方程为 \(x = x(s)\)，其中 \(x(0)=\xi\)。沿着曲线 \(l_\xi\)，\(u\) 的变化情况为

\[ \begin{aligned} \frac{d}{ds}u(x(s),s) = u_x(x(s),x) x'(s) + u_t(x(s),s)\\ = u_x(x(s),s)\, (x'(s) - u(x(s),s)) \end{aligned} \]

若曲线 \(l_\xi:x=x(s)\) 取为如下 ODE 的解

\[ \left\{ \begin{aligned} x'(s) &= u(x(s),s)\\ x(0) &= \xi \end{aligned} \right. \]

那么由上式可知 \(u\) 的值沿曲线不变，这条曲线称为 Burgers 方程的特征线。

由于

\[ x''(s) = \frac{d}{ds} u(x(s),s) = 0 \]

因此 \(l_\xi:x=x(s)\) 为直线，斜率为

\[ x'(s) = x'(s)|_{s=0} = u(x(0),0) = u_0(\xi) \]

最终得到特征线 \(l_\xi\) 的表达式

\[ l_\xi: x(s) = \xi + u_0(\xi) s \]

Newton迭代求解

求解 \(u(x,t)\) ，在不考虑多条特征线相交等复杂情况时，我们只需要找到经过\((x,t)\) 的一条特征线 \(l_\xi\)，然后利用特征线的性质 \(u(x,t) = u_0(\xi)\) 即可。

问题归结为通过如下非线性方程求解得到 \(\xi\) \[ x = \xi + u_0(\xi) t \]

可以使用 Newton 迭代法进行求解：对于一般问题 \(F(\xi) = 0\)，Newton 迭代公式为 \[ \xi^{(k+1)} = \xi^{(k)} - \frac{F(\xi^{(k)})}{F'(\xi^{(k)})} \] 记 \(F(\xi) = \xi -x + u_0(\xi) t\)，Newton 迭代公式为 \[ \xi^{(k+1)} = \xi^{(k)} - \frac{\xi^{(k)}-x+u_0(\xi^{(k)})t}{1 + u_0'(\xi^{(k)})t} \] 选取合适的初值 \(\xi^{(0)}\)，假设 Newton 迭代收敛 \(\lim_{k \to \infty} \xi^{(k)} = \xi^*\)，那么 \[ u(x,t) = u_0(\xi^*) \]

我们还可以根据特征线的性质，直接得到如下关于 \(u\) 的非线性方程 \[ u = u_0(x - u t) \] 求解可以直接获得 \((x,t)\) 位置的数值解 \(u\)

记 \(G(u) = u - u_0(x - u t)\)，为了求解 \(G(u) = 0\)，Newton 迭代公式为 \[ u^{(k+1)} = u^{(k)} - \frac{u^{(k)} - u_0(x - u^{(k)}t)}{1+u_0'(x-u^{(k)}t)t} \] 选取合适的初值 \(u^{(0)}\)，假设 Newton 迭代收敛 \(\lim_{k \to \infty} u^{(k)} = u^*\)，那么 \[ u(x,t) = u^* \]

事实上，这两种做法是等价的，我们只需要考虑如下的变换 \[ \xi^{(k)} = x - u^{(k)} t \]

Newton 迭代求解的关键在于迭代初值的选取，我们可以（并且在激波附近必须）通过特征线的大致趋势或者解的大致形态，选择一个比较好的初值，才能收敛得到正确的结果，具体讨论见下文。

三角函数型初值问题

我们考虑如下初值问题

\[ \left\{ \begin{aligned} u_t + \left(\frac{u^2}2\right)_x &= 0\\ u(x,0) &= \alpha + \beta \sin(wx+\phi) \end{aligned} \right. \]

其中 \(x \in [x_l,x_r]\) ，周期性边界条件，因此系数满足 \(2\pi = w(x_r-x_l)\)，还要求 \(w>0\)，\(\beta \neq 0\)。

只需要进行一次代换，就可以将问题划归为正弦初值问题的求解：\(u(x,t)\) 满足如下形式

\[ u(x,t) = \alpha + \beta \,\widetilde{u}(w x + \phi -\alpha w t,\beta w t) \]

其中 \(\widetilde{u}\) 是在 \([0,2\pi]\) 的正弦初值问题精确解

\[ \left\{ \begin{aligned} \widetilde{u}_t + \left(\frac{\widetilde{u}^2}2\right)_x &= 0\\ \widetilde{u}(x,0) &= \sin(x) \end{aligned} \right. \]

验证一下

\[ \begin{aligned} u &= \alpha + \beta \widetilde{u} \\ u_x &= \beta w \widetilde{u}_x\\ u_t &= -\alpha \beta w \widetilde{u}_x + \beta^2 w \widetilde{u}_t\\ u_t + u u_x &= -\alpha \beta w\widetilde{u}_x + \beta^2 w\widetilde{u}_t + (\alpha + \beta \widetilde{u}) \beta w \widetilde{u}_x \\ &= \beta^2 w(\widetilde{u}_t + \widetilde{u} \widetilde{u}_x ) \end{aligned} \]

对于正弦初值问题，稍微分析一下：

  • 易知精确解 \(\widetilde{u}(x,t)\) 仍然以 \(2\pi\) 为周期，关于 \((\pi,0)\) 中心对称，因此我们只需要关注 \(x\in[0,\pi)\) ，其余区域均可以通过对称性导出。
  • 由特征线的斜率趋势可得，特征线首次交汇的时间

\[ T_b = \frac{-1}{\min u_0'(x)} = 1 \]

  • 精确解 \(\widetilde{u}(x,t)\) 在 \(T_b\) 时刻之后会在 \(x=\pi\) 固定位置产生激波，图像发生间断。根据周期性，在 \(x=(2k+1)\pi\) 位置均有激波产生。

根据代换关系可知，三角函数型初值问题的精确解 \(u(x,t)\) 具有如下性质：

  • 特征线首次交汇的时间为

\[ T_b = \frac{-1}{\min u_0'(x)} = \frac{1}{|\beta w|} \]

  • 激波出现的位置不再固定，而是以速度 \(\alpha\) 移动，具体位置如下（根据周期性可以得到其它激波）

\[ x = \alpha t + \frac{\pi - \phi}w \]

迭代初值选取（近似三角形）

在激波出现之前，由于非线性方程只有一个零点，初值的选取要求其实很宽松，很容易收敛到唯一的零点。但是在激波出现之后，由于此时非线性方程本来就可能有多个零点，而我们希望 Newton 迭代收敛到符合熵解要求的零点，就必须让迭代初值足够得靠近它。

我们首先考虑直接求解 \(u\) 的非线性方程 \[ G(u) = u - u_0(x - u t) \] 我们需要根据解的大致形态来预估一个初值 \(u^{(0)}\)，可以直接将解的形态近似为关于 \((\pi,0)\) 中心对称的，一条边为\(x\)轴的两个直角三角形：

  • 在 \([0,\pi]\) 部分的直角三角形的斜边方程为 \(y = k x\)；
  • 在 \([\pi,2\pi]\) 部分的直角三角形的斜边方程为 \(y=k(x-2\pi)\)。

可以使用如下的初值 \[ u^{(0)}= \left\{ \begin{aligned} & k x, & 0&<x<\pi \\ & k(x-2\pi), & \pi&<x<2\pi \end{aligned} \right. \]

斜率\(k\)的具体选取策略可以有很多，例如：

  • 考虑最值点的移动：初始时刻最值点 \((\pi/2,1)\) 在 \(t\) 时刻会沿着特征线移动到 \((\pi/2+t,1)\)，取它和原点连线的斜率 \(k_1 = 1/(\pi/2+t)\)；
  • 考虑原点处真解的斜率：初始时刻靠近原点的点 \((s,\sin(s))\) 在 \(t\) 时刻会沿着特征线移动到 \((s+\sin(s)t,\sin(s))\)，与原点连线的斜率 \(k(s) = 1/(s/\sin(s)+t)\)，取 \(s \to 0\) 的极限得到 \(k_2 = 1/(1+t)\)；

如下图所示，它们预估的初值都是差不多的，都可以保证收敛。

附上绘图代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
import numpy as np
import matplotlib.pyplot as plt

xleft, xright = 0, 2 * np.pi
num = 100
tend = 2.0

xi_list = np.linspace(xleft, xright, num)
x_list = xi_list + np.sin(xi_list) * tend
u_list = np.sin(xi_list)

fig, axes = plt.subplots()
axes.plot(xi_list, np.sin(xi_list), "b--", alpha=0.3)
axes.plot(x_list, u_list, "b")

axes.plot([np.pi, np.pi], [-1, 1], color="g", linestyle="--")

axes.plot([0, np.pi * 0.5 + 1], [0, 1], color="orange", linestyle="--", lw=0.9)
axes.plot([np.pi * 2, np.pi * 1.5 - 1], [0, -1], color="orange", linestyle="--", lw=0.9)

axes.plot([0, np.pi * 0.5 + tend], [0, 1], color="r", linestyle="--", lw=0.9)
axes.plot([np.pi * 2, np.pi * 1.5 - tend], [0, -1], color="r", linestyle="--", lw=0.9)

axes.plot(np.pi * 0.5 + tend, 1, "rD", markersize=5)
axes.plot(np.pi * 1.5 - tend, -1, "rD", markersize=5)

axes.set_xlim(0, np.pi * 2)
axes.set_ylim(-1.05, 1.05)

axes.spines["right"].set_alpha(0)
axes.spines["top"].set_alpha(0)
axes.spines["bottom"].set_position(("data", 0))
axes.spines["left"].set_position(("data", 0))

axes.set_xticks([0, np.pi, 2 * np.pi])
axes.set_xticklabels(["0", "$\\pi$", "$2\\pi$"])
axes.set_yticks([0, -1, 1])

对于 \(\xi\) 满足的非线性方程的 Newton 迭代求解 \[ F(\xi) = \xi -x + u_0(\xi) t \] 几何含义为找到 \((x,t)\) 所处的，从 \((\xi,0)\) 发出的特征线，根据正弦问题特征线的分布，即使是出现激波之后存在多条特征线交汇， 真实的解所对应的 \((\xi,0)\) 也可以大致进行估计：

  • 对于 \(0<x<\pi\)，必然有 \(0<\xi<\pi\)，即要求的是来自激波左侧的特征线；
  • 对于 \(\pi<x<2\pi\)，必然有 \(\pi<\xi<2\pi\)，即要求的是来自激波右侧的特征线；

对于激波出现之后的迭代可以简单地使用如下的初值 \[ \xi^{(0)}= \left\{ \begin{aligned} & 0, & 0&<x<\pi \\ & 2\pi, & \pi&<x<2\pi \end{aligned} \right. \] 或者基于前面的 \(u^{(0)}\) 的迭代初值讨论，考虑到两者的等价性 \[ \xi^{(0)} = x - u^{(0)} t \] 即 \[ \xi^{(0)}= \left\{ \begin{aligned} & (1-kt) x, & 0&<x<\pi \\ & (1-kt)(x-2\pi)+2\pi, & \pi&<x<2\pi \end{aligned} \right. \] 同样地，可以取 \(k_1 = 1/(\pi/2+t)\) 或 \(k_2 = 1/(1+t)\)。

需要说明的是，上文对迭代初值的讨论都放在了 \([0,2\pi]\) 这个周期，此时绘图和分析更加直观，但是在编程实践中可以考虑 \([-\pi,\pi]\) 这个周期，此时初值的设置就不再需要分段讨论了，可以直接取 \[ u^{(0)} = k x,\quad\text{or}\quad\xi^{(0)} = (1-kt)x \]

迭代初值选取（扫描法）

更一般的做法是扫描法：在一个周期中均匀选取一组初值 \(\{\xi_i\}\)， 然后通过特征线 \[ x = \xi_i + u_0(\xi_i) t \] 将 \(x\!-\!t\) 平面划分为若干区域，直接扫描查找 \(u(x,t)\) 所属的区域，使用附近的特征线所对应的值作为迭代初值。

仍然考虑正弦初值问题，我们基于 \([-\pi,\pi]\) 周期进行扫描。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
import numpy as np
import matplotlib.pyplot as plt

# 在 [-pi,pi] 扫描
def start_point_via_scan(x, t):
    x = np.mod(x + np.pi, 2 * np.pi) - np.pi
    xvec = np.linspace(-np.pi, np.pi, 30)
    uvec = np.sin(xvec)

    condition = (xvec[:, None] + uvec[:, None] * t) >= x
    idx = np.argmax(condition, axis=0)
    return uvec[idx]


xleft = -np.pi
xright = np.pi
xi_list = np.linspace(xleft, xright, 30)
xi_end_list = np.zeros_like(xi_list)

for tend in [0.5, 2.0]:
    fig, axes = plt.subplots(2, 1)
    axes[0].set_xlim(xleft, xright)
    for i, xi in enumerate(xi_list):
        x = np.linspace(xi, xi + np.sin(xi) * tend)
        xi_end_list[i] = x[-1]
        axes[0].plot(x, np.linspace(0, tend), color="b")

    for xi_end in xi_end_list:
        axes[1].axvline(x=xi_end, color="b", linestyle="--")

    x = np.linspace(xleft, xright, 300)
    u = start_point_via_scan(x, tend)
    axes[1].set_xlim(xleft, xright)
    axes[1].plot(x, u, color="r")

    axes[1].plot(x + np.sin(x) * tend, np.sin(x), linestyle="--", color="g")

两个时刻的特征线分布和迭代初值选取如下图所示

由图像可知，这种初值选取策略是 ok 的。

但是这里的周期选择是比较重要的，否则也可能因为特征线交汇和扫描顺序导致错误的结果，例如基于 \([0,2\pi]\) 周期进行扫描。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
# 在 [0,2*pi] 扫描
def start_point_via_scan_2(x, t):
    x = np.mod(x, 2 * np.pi)
    xvec = np.linspace(0, 2 * np.pi, 30)
    uvec = np.sin(xvec)

    condition = (xvec[:, None] + uvec[:, None] * t) >= x
    idx = np.argmax(condition, axis=0)
    return uvec[idx]


xleft = 0
xright = 2 * np.pi
xi_list = np.linspace(xleft, xright, 30)
xi_end_list = np.zeros_like(xi_list)

for tend in [0.5, 2.0]:
    fig, axes = plt.subplots(2, 1)
    axes[0].set_xlim(xleft, xright)
    for i, xi in enumerate(xi_list):
        x = np.linspace(xi, xi + np.sin(xi) * tend)
        xi_end_list[i] = x[-1]
        axes[0].plot(x, np.linspace(0, tend), color="b")

    for xi_end in xi_end_list:
        axes[1].axvline(x=xi_end, color="b", linestyle="--")

    x = np.linspace(xleft, xright, 300)
    u = start_point_via_scan_2(x, tend)
    axes[1].set_xlim(xleft, xright)
    axes[1].plot(x, u, color="r")

    axes[1].plot(x + np.sin(x) * tend, np.sin(x), linestyle="--", color="g")

两个时刻的特征线分布和迭代初值选取如下图所示

由图像可知，这种初值选取策略存在问题。

一般初值问题

考虑如下问题

\[ \left\{ \begin{aligned} u_t + \left(\frac{u^2}2\right)_x &= 0\\ u(x,0) &= u_0(x) \end{aligned} \right. \]

其中 \(x \in [x_l,x_r]\) ，周期性边界条件。

对于一般初值的求解，问题的关键仍然是如下非线性方程的零点求解

\[ x = \xi + u_0(\xi) t \]

但是仅仅通过 Newton 迭代法是无法完美求解的，对于光滑解可以使用扫描法预测一个初值，通常是没有问题的；但是对于间断解则面临着本质上的困难：

  • 间断代表着特征线交汇，非线性方程本来就有两个或多个解，应该选择哪一个解需要基于更多信息进行分析；
  • Newton 迭代对于迭代初值可能很敏感，例如对应在间断附近的两个初值，很可能收敛到差距很远的两个解，或者发散，这在激波附近尤其明显。

Python 求解实现 & 示例

因为三角函数型初值的 Burgers 方程的求解还会频繁使用，因此附上一组 Python 求解函数，包括如下几个函数

  • burgers_sin_exact_solver: 求解函数；
  • burgers_sin_shock: 返回激波信息，包括时间 \(T_b\) 和激波的位置；
  • burgers_sin_show: 调用上述函数求解并绘图，如果存在激波会自动用红线进行标注。

在求解函数中针对非线性方程 \[ G(u) = u - u_0(x - u t) \] 进行Newton迭代，初值的选取策略也很简单 \[ u^{(0)} = k_1 x = \frac{1}{\pi/2+t} x \]

当然也可以将其更换为针对非线性方程 \[ F(\xi) = \xi -x + u_0(\xi) t \] 的Newton迭代，也可以更换其他的初值选取策略。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
def burgers_sin_exact_solver_kernel(x: np.ndarray, t: float, ep: float) -> np.ndarray:
    """
    u(x,t): exact solution of burgers equation when u0(x) = sin(x)
    """

    if t < 0:
        raise ValueError(f"{t=}<0")
    if ep <= 0:
        ep = 1e-6

    # x in [-pi, pi]
    x = np.mod(x + np.pi, 2 * np.pi) - np.pi

    # The choice of initial value ​​in Newton's method is very sensitive!
    k = 1.0 / (np.pi / 2 + t)
    u = k * x

    iter_max: int = 100000
    iter: int = 0
    while iter < iter_max:
        tmp = 1 + np.cos(x - u * t) * t
        du = (u - np.sin(x - u * t)) / tmp
        u = u - du

        iter += 1
        if np.max(np.abs(du)) < ep:
            break

    return u


def burgers_sin_exact_solver(
    x: np.ndarray, t: float, a: float, b: float, w: float, phi: float, ep: float
) -> np.ndarray:
    """
    u(x,t): exact solution of burgers equation when u0(x) = a + b sin(w x + phi)
    """

    return a + b * burgers_sin_exact_solver_kernel(
        w * x + phi - a * w * t, b * w * t, ep
    )


def burgers_sin_shock(
    xleft: float, xright: float, t: float, a: float, b: float, w: float, phi: float
) -> Tuple[float, np.ndarray]:
    """
    Return the earliest time of shock wave appearance and shock positions.

    Tb: Earliest time of shock wave appearance.
    y: Shock positions within the range [xleft, xright].
    """

    tb = np.abs(1 / (b * w))

    if t >= tb:
        y0 = a * t + (np.pi - phi) / w
        tx = 2 * np.pi / w

        y0 = np.mod((y0 - xleft), tx) + xleft
        y = np.arange(y0, xright, tx)
    else:
        y = np.array([])

    return tb, y


def burgers_sin_show(
    xleft: float,
    xright: float,
    xnum: int,
    t: float,
    a: float,
    b: float,
    w: float,
    phi: float,
    ep: float,
    retFigure: bool = False,
):
    x = np.linspace(xleft, xright, xnum)
    u = burgers_sin_exact_solver(x, t, a, b, w, phi, ep)

    fig, axes = plt.subplots()
    tb, y = burgers_sin_shock(xleft, xright, t, a, b, w, phi)

    print(f"u0(x)={a}+{b}*sin({w}*x+{phi})")

    if t >= tb:
        print(f"t({t}) >= tb({tb})")

        if len(y) > 0:
            print(f"Shock wave location: {y}")

            y2 = np.linspace(a - b, a + b, xnum)
            for item in y:
                x2 = item * np.ones(xnum)
                axes.plot(x2, y2, color="r", linestyle="--")
    else:
        print(f"tb = {tb}")

    axes.plot(x, u)

    if retFigure:
        return fig, axes

下面提供几个求解示例。

第一个例子，初值\(u_0(x) = \sin(x)\)，\(t=0.5\)，光滑解情形

1
2
3
4
5
6
7
8
9
10
11
burgers_sin_show(
    xleft=0,
    xright=2 * np.pi,
    xnum=400,
    t=0.5,
    a=0,
    b=1,
    w=1,
    phi=0,
    ep=1e-6,
)

第二个例子，初值\(u_0(x) = \sin(x)\)，\(t=1.3\)，间断解情形

1
2
3
4
5
6
7
8
9
10
11
burgers_sin_show(
    xleft=0,
    xright=2 * np.pi,
    xnum=400,
    t=1.3,
    a=0,
    b=1,
    w=1,
    phi=0,
    ep=1e-6,
)

第三个例子，初值\(u_0(x) = 1 + 2\sin(x)\)，\(t=0.4\)，光滑解情形

1
2
3
4
5
6
7
8
9
10
11
burgers_sin_show(
    xleft=-np.pi,
    xright= np.pi,
    xnum=400,
    t=0.4,
    a=1,
    b=2,
    w=1,
    phi=0,
    ep=1e-6,
)

第四个例子，初值\(u_0(x) = 1 + 2\sin(x)\)，\(t=1.0\)，间断解情形

1
2
3
4
5
6
7
8
9
10
11
burgers_sin_show(
    xleft=-np.pi,
    xright= np.pi,
    xnum=400,
    t=1.0,
    a=1,
    b=2,
    w=1,
    phi=0,
    ep=1e-6,
)

第五个例子，初值\(u_0(x) = 1 + 2\sin(x)\)，\(t=5.0\)，长时间的间断解情形

1
2
3
4
5
6
7
8
9
10
11
burgers_sin_show(
    xleft=-2*np.pi,
    xright= 2*np.pi,
    xnum=400,
    t=5.0,
    a=1,
    b=2,
    w=1,
    phi=0,
    ep=1e-6,
)

第六个例子，周期非\(2\pi\)的情况，初值\(u_0(x) = 0.25 + 0.5\sin(\pi x)\)，\(t=0.7\)，间断解情形

1
2
3
4
5
6
7
8
9
10
11
burgers_sin_show(
    xleft=-1,
    xright=1,
    xnum=400,
    t=0.7,
    a=0.25,
    b=0.5,
    w=np.pi,
    phi=0,
    ep=1e-6,
)

特征线可视化

这里提供基于特征线的可视化函数，可以展示 Burgers 方程的精确解与特征线的关系。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
def burgers_show_line(
    u0: Callable[[np.ndarray], np.ndarray],
    xleft: float,
    xright: float,
    num: int,
    tend: float,
    add_time: bool = True,
):
    xi_list = np.linspace(xleft, xright, num)
    x_list = np.zeros(num)
    u_list = np.zeros(num)

    fig, axes = plt.subplots()

    if add_time:
        axes.plot(xi_list, u0(xi_list), "b", alpha=0.2)

    for i in range(num):
        x = xi_list[i] + u0(xi_list[i]) * tend
        u = u0(xi_list[i])

        xs = np.linspace(xi_list[i], x, num)
        ys = np.linspace(0, tend, num)
        if add_time:
            ys += u

        x_list[i] = x
        u_list[i] = u

        axes.plot(xs, ys, alpha=0.6)

    if add_time:
        axes.plot(x_list, u_list + tend, "r")

使用示例如下，显示初值曲线和当前曲线，以及对应点之间的特征线，尤其注意为了将两个曲线分离并保留特征线，对当前解的曲线的 x 轴相比初值的 x 轴额外提升了\(t\)高度：

  • 初值曲线：\(y=u(x,0)=u_0(x)\)
  • 当前曲线：\(y=u(x,t)+t\)
1
2
3
burgers_show_line(
    u0=lambda s: np.sin(s), xleft=0, xright=2*np.pi, num=100, tend=1.2
)

1
2
3
burgers_show_line(
    u0=lambda s: 1+2*np.sin(s), xleft=-2*np.pi, xright=2*np.pi, num=200, tend=1.0
)

1
2
3
4
5
6
7
8
def gaussian_distribution(x):
    mu = 0
    sigma = 0.6
    return 1 / (sigma * np.sqrt(2 * np.pi)) * np.exp(-0.5 * ((x - mu) / sigma) ** 2)

burgers_show_line(
    u0=gaussian_distribution, xleft=-np.pi, xright=np.pi, num=100, tend=2.0
)

1
2
3
4
5
6
7
8
9
10
11
def jump(x):
    result = np.ones_like(x)
    result[x < 0] = 1
    result[x > 1] = 0
    result[(x >= 0) & (x <= 1)] = 1 - x[(x >= 0) & (x <= 1)]

    return result

burgers_show_line(
    u0=jump, xleft=-2, xright=2, num=100, tend=1.0
)

动态图

顺便提供基于 matplotlib.animation 生成的简单动态图，可以呈现 Burgers 函数解的变化。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
def burgers_show_animate(
    u0: Callable[[np.ndarray], np.ndarray],
    xleft: float,
    xright: float,
    num: int,
    tend: float,
) -> animation.FuncAnimation:
    xi_lists = np.linspace(xleft, xright, num)
    u_lists = u0(xi_lists)
    time_N = 100
    time_dt = tend / time_N

    fig, axes = plt.subplots()

    (line,) = axes.plot(xi_lists, u_lists)
    time_text = axes.text(0.02, 0.95, "", transform=axes.transAxes)

    def update(frame: int):
        x_lists = xi_lists + u0(xi_lists) * frame * time_dt
        line.set_data(x_lists, u_lists)
        time_text.set_text(f"T={frame * time_dt:.2f}")
        return line, time_text

    ani = animation.FuncAnimation(fig=fig, func=update, frames=time_N, interval=20)
    # plt.show()
    return ani

几个不同初值的示例动画如下

1
2
3
4
tmp = burgers_show_animate(
    u0=lambda s: np.sin(s), xleft=0, xright=2 * np.pi, num=100, tend=1.2
)
tmp.save("demo.gif")

1
2
3
4
5
6
7
8
9
def gaussian_distribution(x):
    mu = 0
    sigma = 0.6
    return 1 / (sigma * np.sqrt(2 * np.pi)) * np.exp(-0.5 * ((x - mu) / sigma) ** 2)

tmp = burgers_show_animate(
    u0=gaussian_distribution, xleft=-np.pi, xright=np.pi, num=200, tend=2.0
)
tmp.save("demo.gif")

1
2
3
4
5
6
7
8
9
10
11
12
def jump(x):
    result = np.ones_like(x)
    result[x < 0] = 1
    result[x > 1] = 0
    result[(x >= 0) & (x <= 1)] = 1 - x[(x >= 0) & (x <= 1)]

    return result

tmp = burgers_show_animate(
    u0=jump, xleft=-2, xright=2, num=100, tend=1.0
)
tmp.save("demo.gif")