考察Vlasov-Maxwell系统，证明质量守恒、能量守恒，并探讨粒子方法求解该系统的弱解及能量保持性

考察以下系统：

$\partial _{t}f+v\cdot \nabla _{x}f+E\cdot \nabla _{v}f=0, \quad \partial_{t}E=-\int_{\mathbb{R}^{d}}vf\textrm{d}v$ 。（1）

其中变量有 $t\in\mathbb{R}^{+}$ ， $x\in\mathbb{R}^{d}$ ， $v\in\mathbb{R}^{d}$ ， $f$ 是一个关于 $t,x$ 和 $v$ 的未知函数， $E$ 只关于 $t$ 和 $x$ 。

a)证明质量和能量

$m(t):=\int_{\mathbb{R}^{d}}\int_{\mathbb{R}^{d}}fdvdx$ ，

$e(t):=\frac{1}{2}\int_{\mathbb{R}^{\mathrm{d}}}\int_{\mathbb{R}^{d}}\|v\|^{2}fdxdv+\frac{1}{2}\int_{\mathbb{R}^{d}}|E|^{2}dx$ 。

随时间保持不变。

b）现假设 $E(t,x)$ 已知，则（1）简化为

$\partial_{t}f+v \cdot \nabla_{x} f+E \cdot \nabla_{v} f=0$ 。 (2)

我们将使用粒子方法来求解(2)，初始条件为：

$f(t=0, x, v)=f_{\text {in }}(x, v)$ 。

考虑由 $\left\{x_{p}(t), v_{p}(t)\right\}_{1 \leq p \leq P}$ 表示P个粒子，其中 $x_{p}(t)$ 表示第 $P$ 个粒子的位置， $v_{p}(t)$ 表示其速度，每个粒子的权重为 $\frac{1}{P}$ 。 $x_{p}(t)$ 和 $v_{p}(t)$ 的演化由以下方程决定：

$\frac{d x_{p}}{d t}=v_{p}, \quad \frac{d v_{p}}{d t}=E(x_{p}) 。$

证明

$f^{P}(t, x, v):=\frac{1}{P} \sum_{p=1}^{P} \delta\left(x-x_{p}(t)\right) \delta\left(v-v_{p}(t)\right)$

在分布意义下是(2)的弱解，其中(2)的初始条件是 $f_{0}^{P}=\frac{1}{P} \sum_{p=1}^{P} \delta\left(x-x_{p}(0)\right) \delta\left(v-v_{p}(0)\right)$ 。这里 $\delta$ 是Dirac-Delta函数。

c）回到原始系统（1）。考虑以下粒子方法

$\begin{aligned} & x_{i}^{n+1}=x_{i}^{n}+\Delta t v_{i}^{n} & \\ & v_{p}^{n+1}=v_{p}^{n}-E\left(x_{p}^{n}\right) \Delta t & \\ & E_{i}^{n+1}=E_{i}^{n}-J_{i}^{n} \Delta t \end{aligned}$ 。

这里上标n表示时间步长， $E\left(x_{p}^{n}\right)$ 和 $J_{i}^{n}$ 分别定义如下

$E\left(x_{p}^{n}\right):=\sum_{i} S\left(x_{i}-x_{p}^{n}\right) E_{i}^{n} \Delta x^{d}$ ，

$J_{i}^{n}:=\frac{1}{P} \sum_{p=1}^{P}S\left(x_{i}-x_{p}^{n}\right) v_{p}^{n}$ 。

在这些表达式中， $x_{i}$ 表示具有网格大小 $\Delta x$ 的均匀网格，S是一个满足 $\int_{\mathbb{R}^{d}} S(x) d x=1$ 。这种方法是否保留能量？也就是说，我们是否有

$\frac{1}{2 P} \sum_{p}\left(v_{p}^{n}\right)^{2}+\frac{1}{2} \sum_{i}\left(E_{i}^{n}\right)^{2} \Delta x^{d}=\frac{1}{2 P} \sum_{p}\left(v_{p}^{n+1}\right)^{2}+\frac{1}{2} \sum_{i}\left(E_{i}^{n+1}\right)^{2} \Delta x^{d}$ ？

如果是，证明它。如果不是，构造一个在离散级别上保持能量的离散格式。

证：

a) 我们先证明质量 $m(t)$ 随时间保持不变。

质量 $m(t)$ 定义为

$m(t) = \int_{\mathbb{R}^d} \int_{\mathbb{R}^d} f \, dv \, dx$ 。

对 $m(t)$ 关于时间 $t$ 求导，并使用方程 (1) 中的第一个方程，我们有

$\frac{dm}{dt} = \frac{d}{dt} \int_{\mathbb{R}^d} \int_{\mathbb{R}^d} f \, dv \, dx$ 。

通过分部积分，我们可以将时间导数移到积分内部：

$\frac{dm}{dt} = \int_{\mathbb{R}^d} \int_{\mathbb{R}^d} \partial_t f \, dv \, dx$ 。

根据方程 (1) 中的第一个方程， $\partial_t f + v \cdot \nabla_x f + E \cdot \nabla_v f = 0$ ，因此

$\frac{dm}{dt} = - \int_{\mathbb{R}^d} \int_{\mathbb{R}^d} (v \cdot \nabla_x f + E \cdot \nabla_v f) \, dv \, dx$ 。

使用散度定理，我们可以将 $\nabla_x$ 和 $\nabla_v$ 的积分转换为边界积分。由于 $f$ 在无穷远处趋于零，边界积分为零，因此

$\frac{dm}{dt} = 0$ 。

这表明质量 $m(t)$ 随时间保持不变。

接下来，我们证明能量 $e(t)$ 随时间保持不变。

能量 $e(t)$ 定义为

$e(t) = \frac{1}{2} \int_{\mathbb{R}^d} \int_{\mathbb{R}^d} \|v\|^2 f \, dv \, dx + \frac{1}{2} \int_{\mathbb{R}^d} \|E\|^2 \, dx$ 。

对 $e(t)$ 关于时间 $t$ 求导，并使用方程 (1) 中的两个方程，我们有

$\frac{de}{dt} = \frac{1}{2} \int_{\mathbb{R}^d} \int_{\mathbb{R}^d} 2v \cdot \partial_t v f \, dv \, dx + \frac{1}{2} \int_{\mathbb{R}^d} 2E \cdot \partial_t E \, dx$ 。

根据方程 (1) 中的第一个方程， $\partial_t f = -v \cdot \nabla_x f - E \cdot \nabla_v f$ ，以及第二个方程 $\partial_t E = -\int_{\mathbb{R}^d} v f \, dv$ ，我们可以将 $\partial_t v$ 和 $\partial_t E$ 替换：

$\frac{de}{dt} = \int_{\mathbb{R}^d} \int_{\mathbb{R}^d} (-v \cdot \nabla_x f - E \cdot \nabla_v f) v \, dv \, dx - \int_{\mathbb{R}^d} E \int_{\mathbb{R}^d} v f \, dv \, dx$ 。

使用散度定理和方程 (1) 中的第二个方程，我们可以证明

$\frac{de}{dt} = 0.$

这表明能量 $e(t)$ 随时间保持不变。

b) 现在证明 $f^P(t, x, v)$ 在分布意义下是方程 (2) 的弱解。

首先，我们计算 $f^P(t, x, v)$ 的初始条件：

$f_0^P = \frac{1}{P} \sum_{p=1}^P \delta(x - x_p(0)) \delta(v - v_p(0))$ 。

根据粒子演化方程，我们有

$\frac{d x_p}{dt} = v_p, \quad \frac{d v_p}{dt} = E(x_p)$ 。

对于任意的测试函数 $\phi(t, x, v)$ ，我们计算以下积分：

$\int_{\mathbb{R}^d} \int_{\mathbb{R}^d} \left( \partial_t f^P + v \cdot \nabla_x f^P + E \cdot \nabla_v f^P \right) \phi \, dx \, dv$ 。

利用粒子表示，我们有

$\begin{align*} &\int_{\mathbb{R}^d} \int_{\mathbb{R}^d} \left( \partial_t f^P + v \cdot \nabla_x f^P + E \cdot \nabla_v f^P \right) \phi \, dx \, dv \\ &= \int_{\mathbb{R}^d} \int_{\mathbb{R}^d} \frac{1}{P} \sum_{p=1}^P \left( \partial_t \delta(x - x_p(t)) \delta(v - v_p(t)) + v \cdot \nabla_x \delta(x - x_p(t)) \delta(v - v_p(t)) + E \cdot \nabla_v \delta(x - x_p(t)) \delta(v - v_p(t)) \right) \phi \, dx \, dv \end{align*}$ 。

使用 $\delta$ 函数的性质，我们可以简化上述表达式为

$\begin{align*} &\frac{1}{P} \sum_{p=1}^P \left( \partial_t \phi(t, x_p(t), v_p(t)) + v_p(t) \cdot \nabla_x \phi(t, x_p(t), v_p(t)) + E(x_p(t)) \cdot \nabla_v \phi(t, x_p(t), v_p(t)) \right) \\ &= \frac{1}{P} \sum_{p=1}^P \frac{d}{dt} \phi(t, x_p(t), v_p(t)) \end{align*}$ 。

根据粒子演化方程，我们有

$\begin{align*} \frac{d}{dt} \phi(t, x_p(t), v_p(t)) &= \frac{\partial \phi}{\partial t}(t, x_p(t), v_p(t)) + \frac{d x_p}{dt} \cdot \nabla_x \phi(t, x_p(t), v_p(t)) + \frac{d v_p}{dt} \cdot \nabla_v \phi(t, x_p(t), v_p(t)) \\ &= \frac{\partial \phi}{\partial t}(t, x_p(t), v_p(t)) + v_p(t) \cdot \nabla_x \phi(t, x_p(t), v_p(t)) + E(x_p(t)) \cdot \nabla_v \phi(t, x_p(t), v_p(t)) \end{align*}$ 。

因此，我们有

$\begin{align*} \frac{1}{P} \sum_{p=1}^P \frac{d}{dt} \phi(t, x_p(t), v_p(t)) &= \int_{\mathbb{R}^d} \int_{\mathbb{R}^d} \left( \partial_t f^P + v \cdot \nabla_x f^P + E \cdot \nabla_v f^P \right) \phi \, dx \, dv \\ &= \int_{\mathbb{R}^d} \int_{\mathbb{R}^d} 0 \cdot \phi \, dx \, dv \\ &= 0 \end{align*}$ 。

这表明 $f^P(t, x, v)$ 在分布意义下是方程 (2) 的弱解。

c) 现在考虑原始系统 (1) 的粒子方法。

我们需要验证能量是否在离散时间步长上保持不变。考虑离散化的能量

$E^n = \frac{1}{2P} \sum_{p=1}^P \left(v_p^n\right)^2 + \frac{1}{2} \sum_{i=1}^N \left(E_i^n\right)^2 \Delta x^d,$

其中 $N$ 是网格点的数量。

我们计算 $E^{n+1}$ 并检查是否等于 $E^n$ 。

根据粒子演化方程，我们有

$v_p^{n+1} = v_p^n - E(x_p^n) \Delta t$ 。

因此，

$\begin{align*} \frac{1}{2P} \sum_{p=1}^P \left(v_p^{n+1}\right)^2 &= \frac{1}{2P} \sum_{p=1}^P \left(v_p^n - E(x_p^n) \Delta t\right)^2 &= \frac{1}{2P} \sum_{p=1}^P \left(v_p^n\right)^2 - \frac{1}{P} \sum_{p=1}^P 2v_p^n E(x_p^n) \Delta t + \frac{1}{2P} \sum_{p=1}^P \left(E(x_p^n)\right)^2 (\Delta t)^2 \end{align*}$ 。

接下来，我们考虑电场的更新。假设我们使用显式欧拉方法来更新电场，我们有

$E_i^{n+1} = E_i^n - \Delta t \int_{\mathbb{R}^d} v f^n(x_i, v) \, dv.$

因此，电场能量的变化为

$\begin{align*} \frac{1}{2} \sum_{i=1}^N \left(E_i^{n+1}\right)^2 \Delta x^d &= \frac{1}{2} \sum_{i=1}^N \left(E_i^n - \Delta t \int_{\mathbb{R}^d} v f^n(x_i, v) \, dv\right)^2 \Delta x^d \\ &= \frac{1}{2} \sum_{i=1}^N \left(E_i^n\right)^2 \Delta x^d - \Delta t \sum_{i=1}^N E_i^n \int_{\mathbb{R}^d} v f^n(x_i, v) \, dv \Delta x^d \\ &\quad + \frac{1}{2} \Delta t^2 \sum_{i=1}^N \left(\int_{\mathbb{R}^d} v f^n(x_i, v) \, dv\right)^2 \Delta x^d \end{align*}$ 。

将粒子和电场能量的变化结合起来，我们得到

$\begin{align*} E^{n+1} &= \frac{1}{2P} \sum_{p=1}^P \left(v_p^n\right)^2 - \frac{1}{P} \sum_{p=1}^P 2v_p^n E(x_p^n) \Delta t + \frac{1}{2P} \sum_{p=1}^P \left(E(x_p^n)\right)^2 (\Delta t)^2 \\ &\quad + \frac{1}{2} \sum_{i=1}^N \left(E_i^n\right)^2 \Delta x^d - \Delta t \sum_{i=1}^N E_i^n \int_{\mathbb{R}^d} v f^n(x_i, v) \, dv \Delta x^d \\ &\quad + \frac{1}{2} \Delta t^2 \sum_{i=1}^N \left(\int_{\mathbb{R}^d} v f^n(x_i, v) \, dv\right)^2 \Delta x^d \end{align*}$ 。

为了保持能量，我们需要

$\frac{1}{2P} \sum_{p=1}^P \left(v_p^n\right)^2 - \frac{1}{P} \sum_{p=1}^P 2v_p^n E(x_p^n) \Delta t + \frac{1}{2P} \sum_{p=1}^P \left(E(x_p^n)\right)^2 (\Delta t)^2 + \frac{1}{2} \sum_{i=1}^N \left(E_i^n\right)^2 \Delta x^d - \Delta t \sum_{i=1}^N E_i^n \int_{\mathbb{R}^d} v f^n(x_i, v) \, dv \Delta x^d + \frac{1}{2} \Delta t^2 \sum_{i=1}^N \left(\int_{\mathbb{R}^d} v f^n(x_i, v) \, dv\right)^2 \Delta x^d = E^n$ 。

我们可以看到，如果 $\Delta t$ 足够小，那么 $\left(E(x_p^n)\right)^2 (\Delta t)^2$ 和 $\left(\int_{\mathbb{R}^d} v f^n(x_i, v) \, dv\right)^2 \Delta x^d$ 这两项可以忽略不计，这样能量就会近似保持不变。然而，如果 $\Delta t$ 不是足够小，那么这个方法可能不会保持能量。

为了构造一个在离散级别上保持能量的离散格式，我们可以使用半隐式或隐式时间积分方法。例如，使用隐式欧拉方法更新电场，我们可以得到一个保持能量的离散格式。以下是使用隐式欧拉方法更新电场的步骤：

粒子速度的更新仍然使用显式欧拉方法：

$v_p^{n+1} = v_p^n + \Delta t E(x_p^n)$ 。

电场的更新使用隐式欧拉方法，我们需要解下面的方程来找到 $E_i^{n+1}$ ：

$E_i^{n+1} = E_i^n - \Delta t \int_{\mathbb{R}^d} v f^{n+1}(x_i, v) \, dv$ 。

由于 $f^{n+1}$ 依赖于 $E^{n+1}$ ，我们需要迭代求解 $E^{n+1}$ 和 $f^{n+1}$ 。这通常通过以下步骤完成：

假设 $E^{n+1}$ 的初始猜测，比如 $E^{n+1, 0} = E^n$ 。
使用这个电场值来更新粒子的速度和位置：

$v_p^{n+1, k} = v_p^n + \Delta t E(x_p^n)^k,$

$x_p^{n+1, k} = x_p^n + \Delta t v_p^{n+1, k},$

其中 $k$ 是迭代次数。
使用更新的粒子速度来计算新的 $f^{n+1, k}$ ：

$f^{n+1, k}(x_i, v) = \frac{1}{P} \sum_{p=1}^P \delta(x - x_p^{n+1, k}) \delta(v - v_p^{n+1, k})$ 。
使用新的 $f^{n+1, k}$ 来更新电场：

$E_i^{n+1, k+1} = E_i^n - \Delta t \int_{\mathbb{R}^d} v f^{n+1, k}(x_i, v) \, dv$ 。
重复上述步骤直到 $E^{n+1, k}$ 收敛。

在迭代收敛后，我们得到 $E^{n+1}$ 和 $f^{n+1}$ ，然后可以计算新的能量 $E^{n+1}$ ：

$E^{n+1} = \frac{1}{2P} \sum_{p=1}^P \left(v_p^{n+1}\right)^2 + \frac{1}{2} \sum_{i=1}^N \left(E_i^{n+1}\right)^2 \Delta x^d$ 。

由于隐式方法的稳定性，这种方法通常能够更好地保持能量，即使在较大的时间步长下。这是因为隐式方法考虑了未来时间步的电场值，从而在数值解中引入了稳定性。

要证明这种方法保持能量，我们需要展示在迭代过程中，能量 $E^{n+1}$ 保持不变。这通常涉及到证明迭代过程中的电场和粒子分布函数的更新不会增加总能量。由于这是一个迭代过程，我们可以通过数学归纳法来证明，假设在迭代 $k$ 次后能量保持不变，那么在 $k+1$ 次迭代后能量也将保持不变。这需要详细的数学分析和迭代收敛性的证明，这里不展开。

考察Vlasov-Maxwell系统，证明质量守恒、能量守恒，并探讨粒子方法求解该系统的弱解及能量保持性

推荐阅读更多精彩内容