计算科学 - 在线性 ODE 中利用结构的特例 Runge-Kutta 方法？ - 吾爱随笔录

在线性 ODE 中利用结构的特例 Runge-Kutta 方法？

计算科学颂龙格库塔

2021-12-24 18:49:59

我对以下形式的线性时间相关 ODE 的数值解感兴趣

\dot{y} = A (t) y - R y,

$\dot y = A(t)y - Ry,$

一个好的模型是以下问题 $\mathbb R^2$ ：

A (t) = [\begin{matrix} 0 & - ω (t) \\ ω (t) & 0 \end{matrix}] R = [\begin{matrix} R_{1} & 0 \\ 0 & R_{2} \end{matrix}],

$A(t) = \begin{bmatrix}0 & -\omega(t)\\ \omega(t) & 0\end{bmatrix} \\ R = \begin{bmatrix}R_1 & 0\\ 0 & R_2\end{bmatrix}, \\$

完整的问题在于 $\mathbb R^3$ 并且有点复杂，但我认为上面的示例包含了大部分内容。

该解决方案具有两个主要特征：“振荡”部分和指数衰减。功能 $\omega(t)$ 可以取值高达 $10^6$ , 而衰减因子 $R$ 是在顺序 $10^0$ 到 $10^2$ . 因此，与指数衰减相比，振荡非常迅速。然而， $\omega(t)$ 变化速度很慢，大约为 $10^1$ 到 $10^2$ 每单位时间。

我尝试了经典的 4 阶 Runge-Kutta 方法，效果很好。一个优点是它只需要样本 $A(t)$ 在每个时间步的开始、中间和结束时，这简化了我的实现。

现在我的问题是：是否有可能以某种方式利用该系统的结构来获得更高性能的方法？具体来说，是否可以使用附加结构导出“更好的”龙格-库塔方法（更高阶、更少阶段、更低主误差项常数）？还是使用其他类型的解决方案策略会更好？

编辑：这是我正在使用的方法的 Butcher 画面。为简单起见采用恒定时间步长，因为它将针对许多不同的实例并行运行 $\omega(t)$ 最后在GPU上。我编写了自己的实现，并且验证了错误与时间步长的对数图的斜率为 4。

\begin{array}{ccccc} 0 \\ 1 / 2 & 1 / 2 \\ 1 / 2 & 1 / 2 \\ 1 & 1 \\ 1 / 6 & 1 / 3 & 1 / 3 & 1 / 6 \end{array}

$\begin{array}{c|cccc} 0 & & & & \\ 1/2 & 1/2 & & \\ 1/2& & 1/2 & \\ 1 & & & 1 \\ \hline & 1/6 & 1/3 & 1/3 & 1/6 \end{array}$

1个回答

有多种 RK 方法具有扩展以利用线性。他们都使用某种形式的指数或李群思想（再次指数）来做到这一点。因此，他们通常会采用某种形式的积分因子方法，然后将 Runge-Kutta 方法应用于 IF 变换方程。如果您想查看其中的一些列表，可以查看DifferentialEquations.jl中的方法。

有标准的 RK 方法可以在线性问题上实现超收敛，即达到更高的阶数。DiffEqDevTools.jl实现了我在文献中可以找到的几乎所有表格，并且从中我们发现了以下超收敛案例：

您可以通过在 8,000 行系数表页面中找到它们来挖掘这些方法的论文。看论文，这种超收敛行为似乎是偶然的。但它确实表明可以在线性方程上实现比预期更高的阶数。

现在进入更复杂的方法。在拆分 ODE 求解器部分中，您将找到与以下形式兼容的方法：

u^{'} = A u + f (u, p, t)

$u' = Au + f(u,p,t)$

也称为半线性 ODE。这些被称为指数龙格-库塔方法（指数 RK）。这可能是最常见的用途，但还有其他用途。

如果你只有 $u' = Au$ 那么当然 $u(t) = \exp(TA)u(0)$ 但是您不希望只使用大多数语言中内置的缩放和平方方法，因此您希望使用利用 Krylov 子空间进行大指数时间向量计算的东西，因此该LinearExponential方法。

另一个常见的情况是专门研究的方法 $u' = A(t)u$ ，这是大多数 Magnus 方法所擅长的。您可以找到该类别的方法列表作为状态无关求解器。您写下的问题属于 Magnus 方法，因为您可以 $\hat{A}(t) = A(t) - R$ .

然后是李群方法，如龙格-库塔-蒙特-卡斯方法，专门研究形式 $u' = A(u)u$ . 那么当然有一些方法允许时间和状态依赖，即 $u' = A(u,t)u$ ，特别是一些 Crouch-Grossman Runge Kutta 方法，但随着您的专业化程度较低，它往往会获得越来越少的好处。

与任何DifferentialEquations.jl 方法一样，?alg如果您想了解更多信息，您可以从Julia REPL 中对算法进行挖掘以挖掘原始来源。

我想留下最后一个细节，尽管人们倾向于问这种事情。假设您有一个仿射方程：

u^{'} = A (t) u + b (t)

$u' = A(t)u + b(t)$

你还能用 Magnus 方法利用线性结构吗？答案是肯定的。您需要做的是利用仿射变换只是更大空间中的线性变换这一事实。即你定义 $\hat{u}$ 这样 $\hat{u} = [u;1]$ , 所以你将 1 附加到你的原始空间，然后 $\hat{A}(t)$ 有一个块是 $A$ ，在底部和列上附加一行 0 $b$ 附在右边。如果你计算出来，你会看到方程的扩展部分有一个初始条件 1 和 0 导数，所以 $\hat{A}\hat{u} = A(t)u + b(t)$ 在非常数部分。这意味着您还可以将仿射变换编码为一种形式，以使用这些专门的 RK 方法。

所以 Magnus、Crouch-Grossman 类型的 RK、Runge-Kutta-Munthe-Kaas、指数 Runge-Kutta、Cayley 方法等等。有很多方法可以利用这个属性，并且有一些库可以实现它们。特别是，DifferentialEquations.jl 使用专门的矩阵指数方法，与 GPU 兼容，自动构造伴随等。也就是说，这一切还远未完成，还有很多事情可以做。但这是我所知道的最完整的清单。有关这些方法的背景的更多信息，您可能需要查看Hairer 的几何数值积分一书，该书对此处描述的方法的子集进行了很好的讨论。

其它你可能感兴趣的问题

上一篇Lustre 文件系统的性能下一篇为矩阵的低秩更新廉价地重新计算特征值和特征向量