计算科学 - 刚性方程求解器如何工作？ - 吾爱随笔录

刚性方程求解器如何工作？

计算科学 pde 有限差分数字刚性

2021-12-06 17:00:09

我试图了解如何解决刚性微分方程。

例如方程，

\frac{\partial y}{\partial t} = α \frac{\partial^{2} y}{\partial z^{2}}

$\frac{\partial y}{\partial t} = \alpha\frac{\partial ^2 y}{\partial z^2}$

可以使用 ode 求解器通过离散化第 i 个节点的 z 方向来求解

\frac{d y_{i}}{d t} = α \frac{y_{i + 1} - 2 y_{i} + y_{i - 1}}{Δ z^{2}}

$\frac{dy_i}{dt} = \alpha\frac{y_{i+1}-2y_i + y_{i-1}}{\Delta z^2}$

这个 ode 通常使用 ode15s 来解决。据我了解，刚性求解器使用数值微分形式，如反向微分公式（BDF）来近似导数。考虑一个两步 BDF 公式，

\dot{y} (t_{n}) = \frac{\frac{3}{2} \cdot y_{n} + (- 2) \cdot y_{n - 1} + \frac{1}{2} y_{n - 2}}{Δ t}

$\dot y(t_n)= \frac{\frac{3}{2}\cdot y_n + (-2)\cdot y_{n-1} + \frac{1}{2}y_{n-2}}{\Delta t }$ .

我不知道如何从这里开始。我没有使用求解器，而是尝试实现上述 BDF 以了解求解器的工作原理。

我想就如何从这里开始寻求建议。

编辑：等式（2）和（3）等式

\frac{\frac{3}{2} \cdot y_{i}^{n} + (- 2) \cdot y_{i}^{n - 1} + \frac{1}{2} y_{i}^{n - 2}}{Δ t} = α \frac{y_{i + 1}^{n} - 2 y_{i}^{n} + y_{i - 1}^{n}}{Δ z^{2}}

$\frac{\frac{3}{2}\cdot y^n_i + (-2)\cdot y^{n-1}_i + \frac{1}{2}y^{n-2}_i}{\Delta t }= \alpha\frac{y^n_{i+1}-2y^n_i + y^n_{i-1}}{\Delta z^2}$ 这里，n 代表时间，i 代表空间。

3个回答

您在编辑中设置方程式的方式允许您使用 thomas 算法求解。如果您查找 TDMA 的算法，您会看到该等式将创建一个包含 (n) 时间步的系数的三对角矩阵 A，以及根据您之前的时间步 (n-1, n-2) 的右侧 b解值。这将允许您轻松解决 BDF1 时间步长问题。托马斯算法的维基页面在这里： https ://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm

如果您想解释 BDF1 的工作原理，您可以这样想，我将针对稳态欧拉 CFD 代码进行解释。在稳态欧拉中，我们有保守变量的向量， $u$ ，我们的残差算子对应于离散和集成的 PDE，它必须为零。

R (u) = 0

$R(u) = 0$

我们从 $u_n$ , 为此 $R(u_n)$ 不等于零，并且一个 $R(u_{n+1})$ 我们希望它等于零。如果我们采用泰勒级数展开 $R(u_{n+1})$ 关于 $u_n$ 我们得到：

R (u_{n + 1}) = R (u_{n}) + \frac{\partial R}{\partial u_{n}} Δ u = 0

$R(u_{n+1}) = R(u_n) + \frac{\partial R}{\partial u_n}\Delta u = 0$ 然后我们可以将左项带入 RHS 并得到：

\frac{\partial R}{\partial u_{n}} Δ u = - R (u_{n})

$\frac{\partial R}{\partial u_n}\Delta u = -R(u_n)$ 现在，通过求解这个线性系统，我们实际上不会一步一步求解非线性控制方程，因为线性泰勒级数展开不足以快速求解非线性方程，但我们可以多次重复此迭代以求解稳态状态问题，您可能会注意到这与牛顿方法的相似之处更多。希望这有帮助！

最新的时间应该是未知的，因此可以通过将 lhs 移动到 rhs 来将隐式方法重塑为求根方程组。

不过，您似乎将 y 离散化是错误的，空间离散化需要与时间离散化分开进行，因此如果您保留整个历史，答案是矩阵。（注意：您只显示了时间导数的表达式，而不是积分或泰勒展开式的表达式，并且空间模板的 lhs 不应该是空间导数吗？）

使用ode15s

[t,y] = ode15s(@(t,z) fun(t,z,alpha,dz),tspan,y0);

和

function dy = fun(t,y,alpha,dz)
dy = zeros(size(y));
dy(1) =  % b.c.
dy(2:end-1) = alpha*(y(3:end) -2*y(2:end-1) + y(1:end-2) )/dz^2
dy(end) = % b.c.

哪里dy(1)和dy(end)取决于边界条件。

如果您想实现自己的数值方法，您可以搜索任何方法来求解线性 ODE 系统。

其它你可能感兴趣的问题

上一篇机器学习中具有特殊条件的激活函数下一篇Runge-Kutta，所有节点为 n+1 或零权重