Burgers方程精确解
一维 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)\)的零点，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 | import numpy as np |
对于 \(\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,\,\text{or}\,\xi^{(0)} = (1-kt)x \]
一般初值问题
考虑如下问题
\[ \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迭代对于迭代初值可能很敏感，例如对应在间断附近的两个初值，很可能收敛到差距很远的两个解，或者发散，这在激波附近尤其明显。
一个可选的做法是先采用Godunov格式计算一个数值解，然后将其作为Newton迭代的初值，进一步获取精确解。
Python 求解实现 & 示例
因为三角函数型初值的 Burgers 方程的求解还会频繁使用，因此附上一组 Python 求解函数，包括如下几个函数
burgers_sin_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 | import numpy as np |
下面提供几个求解示例。
第一个例子，初值\(u_0(x) = \sin(x)\)，\(t=0.5\)，光滑解情形
1 | burgers_sin_show( |
第二个例子，初值\(u_0(x) = \sin(x)\)，\(t=1.3\)，间断解情形
1 | burgers_sin_show( |
第三个例子，初值\(u_0(x) = 1 + 2\sin(x)\)，\(t=0.4\)，光滑解情形
1 | burgers_sin_show( |
第四个例子，初值\(u_0(x) = 1 + 2\sin(x)\)，\(t=1.0\)，间断解情形
1 | burgers_sin_show( |
第五个例子，初值\(u_0(x) = 1 + 2\sin(x)\)，\(t=5.0\)，长时间的间断解情形
1 | burgers_sin_show( |
第六个例子，周期非\(2\pi\)的情况，初值\(u_0(x) = 0.25 + 0.5\sin(\pi x)\)，\(t=0.7\)，间断解情形
1 | burgers_sin_show( |
特征线可视化
这里提供基于特征线的可视化函数，可以展示 Burgers 方程的精确解与特征线的关系。
1 | def burgers_show_line(u0, xleft, xright, num, tend): |
使用示例如下，显示初值曲线和当前曲线，以及对应点之间的特征线，尤其注意为了将两个曲线分离并保留特征线，对当前解的曲线的 x 轴相比初值的 x 轴额外提升了\(t\)高度：
- 初值曲线：\(y=u(x,0)=u_0(x)\)
- 当前曲线：\(y=u(x,t)+t\)
1 | burgers_show_line( |
1 | burgers_show_line( |
1 | def gaussian_distribution(x): |
1 | def jump(x): |
动态图
顺便提供基于 matplotlib.animation 生成的简单动态图，可以呈现 Burgers 函数解的变化
1 | import numpy as np |
几个不同初值的示例动画如下
1 | tmp = burgers_show_animate( |
1 | def gaussian_distribution(x): |
1 | def jump(x): |
补充
MATLAB 求解实现
下面是基于MATLAB的简化实现 1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23% Exact solution of Burgers' equation with 2pi periodic boundary condition
% u_t + u*u_x = 0, x in [-pi, pi]
function u = BurgersExactSolution(x,t,alpha,beta)
% u(x,0) = alpha + beta*sin(x)
u = alpha + beta * BugersSinNewton(x - alpha*t, beta*t);
end
function u = BugersSinNewton(x, t, ep)
% u(x,0) = sin(x)
x = mod(x + pi, 2*pi) - pi; % x in [-pi,pi]
iter = 1;
u = x/(pi/2+t); % u0
while iter < 1e5
du = (u-sin(x-u*t))./(1 + cos(x - u*t)*t);
u = u - du;
if max(abs(du)) < 1e-10
break
end
iter = iter + 1;
end
end
Cpp 求解实现
同时附上 C++
版本的三角函数型初值问题的求解实现，封装为一个类BurgersExact，包括：
eval：求解函数；initial_value：获取初值；get_tb：获取激波出现的最早时刻 tb；
没有提供作图或数据导出功能。
1 |
|