Hamilton–Jacobi 方程特征线法求精确解

与 Burgers 方程为代表的守恒律方程类似，Hamilton–Jacobi 方程也可以利用特征线方法进行精确求解。

一维 Hamilton–Jacobi 方程如下

\[ \left\{ \begin{aligned} u_t + H\left(u_x\right) &= 0\\ u(x,0) &= u_0(x) \end{aligned} \right. \]

其中定义域 \(x \in \mathbb{R}\)，通常假定初值有周期性，精确解此时也具有周期性。为了简化记号，引入符号 \(p = u_x\)。

我们将考虑两种具体情况：

• \(H(p) = \frac12 (p+\alpha)^2\)，其中 \(\alpha > 0\) 为常数，此时 \(H(p)\) 是一个严格凸函数；
• \(H(p) = - \cos(p + \alpha)\)，此时 \(H(p)\) 是一个非凸非凹的函数。

特征线方程

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

\[ \begin{aligned} \frac{d}{ds}u_x(x(s),s) = u_{xx}(x(s),s) x'(s) + u_{xt}(x(s),s)\\ \end{aligned} \]

对方程求导可得

\[ u_{tx} + H'(u_x) u_{xx} = 0 \]

实际上，如果记 \(p=u_x\)，可知 \(p\) 满足一个标准的双曲守恒律方程 \(p_t + H(p)_x = 0\)。

代入得到

\[ \begin{aligned} \frac{d}{ds}u_x(x(s),s) = u_{xx}(x(s),s) \Big[ x'(s) - H'(u_x(x(s),s))\Big]\\ \end{aligned} \]

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

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

那么由上式可知 \(u_x\) 的值沿曲线不变，这条曲线称为 H-J 方程的特征线，并且 \(l_\xi:x=x(s)\) 为直线，斜率为

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

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

\[ l_\xi: x(s) = \xi + H'(u_x(\xi,0)) s \]

考虑 \(u\) 沿着特征线 \(l_\xi\) 的变化

\[ \begin{aligned} \frac{d}{ds}u(x(s),s) ={}& u_{x}(x(s),s) x'(s) + u_{t}(x(s),s)\\ ={}& u_{x}(x(s),s) H'(u_x(x(s),s)) - H(u_x(x(s),s)) \\ ={}& u_{x}(x(0),0) H'(u_x(x(0),0)) - H(u_x(x(0),0)) \\ ={}& u_{x}(\xi,0) H'(u_x(\xi,0)) - H(u_x(\xi,0)) \end{aligned} \]

因此 \(u\) 沿着特征线的变化是均匀的，可以表示为

\[ u(x(s),s) = u(\xi,0) + \Big[ u_{x}(\xi,0) H'(u_x(\xi,0)) - H(u_x(\xi,0)) \Big] s \]

可以引入 \(H(p)\) 的对偶函数 \(L(p) = p H'(p) − H(p)\)，\(u\) 可以表示为

\[ u(x(s), s) = u(\xi,0) + L(u_x(\xi,0)) s \]

整理一下：

• 从初始位置 \(x(0) = \xi\) 发出的特征线 \(l_\xi\) 为一条直线

\[ l_\xi: x(s) = \xi + H'(u_x(\xi,0)) s \]

• 沿着特征线 \(l_\xi\)，\(p = u_x\) 保持不变

\[ p(x(s),s) = u_x(x(s),s) = u_x(\xi,0) \]

• 沿着特征线 \(l_\xi\)，\(u\) 均匀变化

\[ \begin{aligned} u(x(s),s) ={}& u(\xi,0) + \Big[ u_{x}(\xi,0) H'(u_x(\xi,0)) - H(u_x(\xi,0)) \Big] s \\ ={}& u(\xi,0) + L(u_x(\xi,0)) s,\quad (L(p) := p H'(p) − H(p)) \end{aligned} \]

与守恒律方程不同，H-J 方程的特征线交汇会导致 \(p = u_x\) 产生多值， 此时 \(u_x\) 不连续，\(u\) 的函数曲线会出现尖点。

对于特征线的交汇和尖点的位置追踪，直接从 \(p\) 满足的守恒律方程入手会更加简单

\[ p_t + H(p)_x = 0, \quad p(x,0) = p_0(x) \]

对应的特征线方程为 \[ l_\xi: x(s) = \xi + H'(p_0(\xi)) s \]

分析易得，初值的光滑性可以保证特征线在短时间内不会交汇，发生交汇的临界时刻为 \[ T_b = \frac{-1}{\min_x H''(p_0(x)) p_0'(x,0)} \]

Newton迭代求解

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

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

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

提供一份求解代码，需要提供迭代初值函数

import numpy as np
import matplotlib.pyplot as plt
from typing import Callable

def HJ_exact(
    x: np.ndarray,
    t: float,
    u0: Callable[[np.ndarray], np.ndarray],
    u0_x: Callable[[np.ndarray], np.ndarray],
    u0_xx: Callable[[np.ndarray], np.ndarray],
    Hp: Callable,
    dHdp: Callable,
    dH2d2p: Callable,
    init: Callable,
    ep: float,
) -> np.ndarray:
    """
    u(x,t): exact solution of the H-J equation
    """

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

    def L(p):
        return p * dHdp(p) - Hp(p)

    xi = init(x, t)

    iter_max: int = 100000
    iter: int = 0
    while iter < iter_max:
        tmp = 1 + dH2d2p(u0_x(xi)) * u0_xx(xi) * t
        dxi = (xi - x + dHdp(u0_x(xi)) * t) / tmp
        xi = xi - dxi

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

    return u0(xi) + L(u0_x(xi)) * t

例一

取计算区域 \([-1,1]\)，初值 \[ u_0(x) = -\cos(\pi x), \quad p_0(x) = \pi \sin(\pi x). \]

考虑 \(H(p) = \frac12(p+\alpha)^2\)（\(\alpha > 0\)），计算可得 \[ H''(p) =1,\quad p_0'(x) = \pi^2 \cos(\pi x) \]

因此临界时刻的精确值为

\[ T_b = \frac{-1}{ \pi^2 \min_x \cos(\pi x)} = \frac{1}{ \pi^2} \approx 0.101 \]

取 \(q = p + \alpha\)，可以得到 \(q\) 满足的守恒律方程

\[ q_t + G(q)_x = 0, \quad q(x,0) = \alpha + p_0(x) = \alpha + \pi \sin(\pi x). \]

其中 \(G(q) = H(q-\alpha) = \frac12 q^2\) 为标准的 Burgers 方程。

由 Burgers 方程的相关结论可知，\(q\) 的间断位置（同样也是 \(p = u_x\) 的间断位置）

\[ x = \alpha t + 1 \]

根据周期性，所有的间断位置为 \(x = \alpha t + 2k + 1\)（\(k \in \mathbb{Z}\)）。

对特征线的分布趋势进行绘图分析可知，整个时空定义域可以被间断位置所对应的直线 \(x_k(t) = \alpha t + 2k+1\) 划分为不相交区域，每个区域内的特征线都向两侧吸引。 特征线分布趋势决定了我们在 Newton 迭代时应当如何选取初值，一个可行的策略为：对于 \((x,t)\)，取 \(\xi^{(0)}\) 为 \(x-\alpha t\) 最接近的偶数。

求解代码如下

def nearest_even(x):
    x = np.asarray(x)
    x_rounded = np.round(x)
    odd_mask = x_rounded % 2 != 0
    lower_even = x_rounded - 1
    upper_even = x_rounded + 1
    nearest = np.where(
        odd_mask,
        np.where(
            np.abs(x - lower_even) <= np.abs(x - upper_even), lower_even, upper_even
        ),
        x_rounded,
    )
    return nearest if isinstance(x, np.ndarray) else nearest.item()


x = np.linspace(-1, 1, 200)
u = HJ_exact(
    x=x,
    t=1.5 / np.pi**2,
    u0=lambda s: -np.cos(np.pi * s),
    u0_x=lambda s: np.pi * np.sin(np.pi * s),
    u0_xx=lambda s: np.pi**2 * np.cos(np.pi * s),
    Hp=lambda p: ((p + 1) ** 2) / 2,
    dHdp=lambda p: p + 1,
    dH2d2p=lambda p: 1,
    init=lambda x, t: nearest_even(x - t),
    ep=1e-8,
)

plt.plot(x, u)

例二

取计算区域 \([-1,1]\)，初值 \[ u_0(x) = -\cos(\pi x), \quad p_0(x) = \pi \sin(\pi x). \]

考虑 \(H(p) = - \cos(p + \alpha)\)（\(\alpha > 0\)），计算可得 \[ H''(p) = \cos(p + \alpha) = \cos(\pi \sin(\pi x) + \alpha),\quad p_0'(x) = \pi^2 \cos(\pi x) \]

因此临界时刻为

\[ T_b = \frac{-1}{\pi^2 \min_x \cos(\pi \sin(\pi x) + \alpha)\cos(\pi x)} \]

并没有简单的显式表达，对于 \(\alpha = 1\) 的情况，数值求解可得 \(T_b \approx 0.106\)。

取 \(q = p + \alpha\)，可以得到 \(q\) 满足的守恒律方程

\[ q_t + G(q)_x = 0, \quad q(x,0) = \alpha + p_0(x) = \alpha + \pi \sin(\pi x). \]

其中 \(G(q) = H(q-\alpha) = - \cos(q)\)。进一步取 \(w = \sin(q)\)，可以得到

\[ w_t + w w_x = 0, \quad w(x,0) = \sin(\alpha + \pi \sin(\pi x)). \]

满足标准的 Burgers 方程。

无论是否代换，对于特征线的分布趋势，我们都无法给出精确的描述，只有一些短时间内的粗糙估计：

• 特征线的走向取决于 \(w_0(\xi)\) 的正负，比较混乱。
• 在 \(t \to T_b^-\)，特征线会在 \(x_1=\xi_1 + w_0(\xi_1) T_b\) 位置汇聚，产生第一个激波；
• 在 \(t \to T_2^-\)，特征线会在 \(x_2=\xi_2 + w_0(\xi_2) T_2\) 位置汇聚，产生第二个激波。

其中 \[ \begin{gathered} T_b \approx 0.106285, \quad \xi_1 \approx -0.907423, \quad x_1 \approx -0.896904 \\ T_2 \approx 0.13088, \quad \xi_2 \approx 0.201604, \quad x_2 \approx 0.238053 \end{gathered} \]

并且对于激波的移动速度似乎难以确定（等面积法？），因此并没有什么好的初值选取方案，选择使用通用的扫描法：预先均匀采样计算一组特征线，然后挑选与 \((x,t)\) 最近的一个作为初值即可。（但是实验发现在尖点附近的效果可能不好，还是因为迭代初值选取得不够好）

求解代码如下

def start_point_via_scan(x, t):
    xstar = 0
    x = np.mod(x - xstar, 2) + xstar
    xvec = np.linspace(xstar, xstar + 2, 40)
    uvec = np.sin(1 + np.pi * np.sin(np.pi * xvec))

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


x = np.linspace(-1, 1, 200)
u = HJ_exact(
    x=x,
    t=1.5 / np.pi**2,
    u0=lambda s: -np.cos(np.pi * s),
    u0_x=lambda s: np.pi * np.sin(np.pi * s),
    u0_xx=lambda s: np.pi**2 * np.cos(np.pi * s),
    Hp=lambda p: -np.cos(p + 1),
    dHdp=lambda p: np.sin(p + 1),
    dH2d2p=lambda p: np.cos(p + 1),
    init=lambda x, t: start_point_via_scan(x, t),
    ep=1e-9,
)

plt.plot(x, u)