计算科学 - 在不连续性上未能与 DotNumerics 中的 radau5 实现集成 - 吾爱随笔录

在不连续性上未能与 DotNumerics 中的 radau5 实现集成

计算科学数值分析颂刚性

2021-12-26 00:10:19

概括

我们正在尝试解决一些物质排放到空气中的模型。在这些模型中，排放在某个点停止。

我们正在使用radau5的DotNumerics实现

我们在寻找正确解决方案时遇到了问题。如果排放停止，空气浓度应该会下降，但情况并非总是如此。

我们能够用一个简单的模型重现这个问题，我们也可以找到分析解决方案。此外，我们检查了 radau5 源代码，发现了一个看起来像优化的测试。由于这个测试，雅可比的许多评估被跳过。如果我们去掉这个测试，结果会好很多。

我们还使用Press et all., Numerical Recipes 2nd edition的隐式四阶 Runge Kutta 来解决我们最初的问题，这种方法没有显示这个问题。我们没有时间对这个问题进行测试。

我们想知道：

这是什么测试？有相关文件吗？
我们可以安全地删除测试以跳过对雅可比的评估吗？
我们能否期望 radau5 能够应对不连续性，或者我们是否应该分别整合这两个轨迹？

详细资料

该模型由两个隔间组成。最初，隔间 I 包含排放到隔间 II 的物质。通过通风从隔间 II 中取出产品。

$$ \begin{array}{rclr} \frac{dC_I}{dt}&=&-k_{12}(C_I- C_{II}) \\ \frac{dC_{II}}{dt}&=&k_ {12}(C_I- C_{II})-k_{02}C_{II} & t < T \\ \frac{dC_{II}}{dt}&=&-k_{02}C_{II} & t > T \\ C_I (0)&=&1 \\ C_{II}(0)&=&0 \end{数组} $$

我们找到了这些分析解决方案：

$$ \begin{array}{rclr} C_{II}(t)&=&C_{II}(T)e^{-k_{02}(tT)} \\ &=&\frac{C_0}{\ sqrt{D}}k_{12}(e^{λ_+ T}-e^{λ_-T})e^{-k_{02}(tT)} & t < T \\ C_{II}(t )&=&\frac{C_0}{\sqrt{D}} k_{12} (e^{λ_+t}-e^{λ_- t}) & t > T \\ D&=&4k_{12}^ 2+k_{02}^2 \\ λ_±&=&-k_{12}-\frac{k_{02}}{2}±\frac{\sqrt{D}}{2} \end{数组} $$

我们在 C# 中实现了这个模型。

公差为 0.01 的数值解如下：

原始 radau5 和相对容差的 jsfiddle：0.01

在这种情况下，Jacobian 仅被评估两次。

翻阅radau5代码，我们发现了一个测试：

if (THETA <= THET) goto LABEL20;

在 .\DotNumerics\ODE\Radau5\radau5.cs 的第 1663 行，

IF (THETA.LE.THET) GOTO 20

在http://www.unige.ch/~hairer/prog/stiff/radau5.f的第 1073 行。

如果我们应用修复，我们会得到：

原始 radau5 和相对容差的 jsfiddle：0.01

将容差降低到 0.001 会产生更好的结果，无论有无修正。

用于原始 radau5 和相对容差的 jsfiddle：0.001 用于修改的 radau5 和相对容差的 jsfiddle：0.001

然而，如果没有修正，雅可比仍然只被评估了 3 次，所以这可能是一个巧合。

在 GitHub 上提供了一个演示这一点的c# 解决方案。

2个回答

首先，我无法判断 Radau5 的 DotNumerics 实现的质量，我只能假设它是 Hairer 和 Wanner 的 Fortran 原始版本的直接端口。他们的 Radau5 实现的完整描述在他们关于求解刚性微分方程的书中。

您所指的测试是检查是否应重新计算雅可比行列式。如果您不提供显式雅可比行列式，Radau 将尝试通过有限差分来估计此雅可比行列式。Radau 代码基于牛顿迭代计算参数“theta”。如果该值大于某个阈值（“THET”），则重新计算雅可比。用户可以选择通过 WORK(3) 变量设置 THET 的值（至少在 Fortran 实现中）。所以不要更改代码中的任何内容！使用此选项可以找到合适的 THET 值。根据 Fortran 实现，WORK(3) 的默认值为 0.001。您可以检查在 DotNumerics 实现中是否保留了此默认值。同样，在 FORTRAN 实现中，

其次，我按照您在 Python 中使用 Assimulo 的描述运行了您的问题。请参阅下面的代码。解决方案变得更粗糙，但我没有观察到跳跃

第三，如果我理解正确，您明确地知道您的“交叉”点 $T$，即它不依赖于问题/模拟。那么，为什么要通过在 ODE 的右侧出现不连续性来让 Radau 的生活变得艰难呢？为什么不简单地积分到$t = T$，然后将这些值用作下一个问题的初始值？

最后，如果你希望有“问题相关”的不连续性，我建议使用 CVODE 代码，它实现了一种非常好的方式来处理状态变化。例如，请参阅我对最近关于如何在满足特定条件时停止模拟的问题的回答以及 Assimulo 网站上描述摆锤撞墙的摆锤问题的示例。

代码示例

import numpy as np
from assimulo.problem import Explicit_Problem
from assimulo.solvers import Radau5ODE
import matplotlib.pyplot as plt
%matplotlib inline

C0 = 1.0
T = 20.0
tf = 60.0
k02 = 0.005
k12 = 0.1

def rhs(t,y):

    ydot = np.zeros_like(y)

    if t <= T:
        ydot[0] = k12 * (y[1] - y[0]) - k02 * y[0]
        ydot[1] = -k12 * (y[1] - y[0])
    else:
        ydot[0] = -k02 * y[0]
        ydot[1] = 0.0

    return ydot

y0 = np.array([C0,0.0])

model = Explicit_Problem(rhs,y0)
simulation = Radau5ODE(model)
simulation.rtol = 1e-2

t, y = simulation.simulate(60.)
simulation.reset()

运行它会给出以下输出

    Final Run Statistics: --- 

 Number of steps                                 : 9
 Number of function evaluations                  : 48
 Number of Jacobian evaluations                  : 2
 Number of function eval. due to Jacobian eval.  : 4
 Number of error test failures                   : 0
 Number of LU decompositions                     : 9

Solver options:

 Solver                  : Radau5 (explicit)
 Tolerances (absolute)   : 0.01
 Tolerances (relative)   : 0.01

Simulation interval    : 0.0 - 60.0 seconds.
Elapsed simulation time: 0.0007725144052983524 seconds.

并绘制两个浓度给出

为了完整起见，这是请求的绝对和相对容差 1e-12 的图。

这与您观察到的不同：在此模拟中没有减少（我从您的 Github 代码中获取了系数的值）。

我们能否期望 radau5 能够应对不连续性，或者我们是否应该分别整合这两个轨迹？

微分方程方法不能轻易处理这样的不连续性。如果你跨过一个不连续点，你不能证明你不会有订单损失。事实上，您通常会损失准确性。因此，您要确保“准确”地击中每个不连续点，即直接踩到它（或将两个部分分别整合）。

对此的一个案例研究是延迟微分方程求解器，其中传播了初始不连续性，延迟方程求解器的目标是使用带有插值的 ODE 方法，以您跟踪和命中不连续性的特定方式。由此可以很好地量化不跟踪不连续性的错误成本。

其它你可能感兴趣的问题

上一篇从 FEM 软件中提取系统矩阵下一篇我可以使用多重网格来解决不是由微分方程产生的线性代数问题吗？