计算科学 - Euler正演方法的步长和稳定性 - 吾爱随笔录

Euler正演方法的步长和稳定性

计算科学稳定显式方法

2021-12-26 06:26:06

我正在尝试使用欧拉前向方法计算为以下非线性 IVP 提供稳定性的最大步长：

$u'(t) = -200tu(t)^2,\qquad u_0 = 1, \qquad t\in [0,3]$ ,

解析解为 $u = (1+100t^2)^{-1}$ （见下图）。从线性稳定性分析，即通过求解 $u' = \lambda u, \qquad u_0 = 1, \qquad \lambda <0$ ，可以证明欧拉前锋的稳定区域是 $|1+\lambda h|\le1$ ，在哪里 $h$ 是步长。所以步长必须是 $h<2/|\lambda|$ .
我意识到这仅适用于线性问题，但是在解决一般 IVP（包括非线性项）时：

$u'(t) = f(t,u(t)), \qquad u(t_0)=u_0$

计算特征值 $\lambda_i(t)$ 的 $A=f_{u_i}(t,u(t))$ ，在哪里 $f_{u_i}$ 是雅可比矩阵，当使用所有的最小值时，也应该导致合理的步长 $\lambda_i$ . 也就是说，如果 $u(t)$ 是渐近稳定的。
所以下一步是计算雅可比矩阵（在这种情况下只有一个条目）：
$f_{u_i}=-\frac{400t}{1+100t^2}=\lambda$ . 下确界 $\lambda$ 在给定的间隔上是 $\lambda=-20$ 因此 $h = 2/20$ . 下图显示 $\lambda(t)$ . 然而，随着模型的爆炸，步长似乎太大了。稳定似乎在某个地方

$h=2/29$ . 所以我有两个问题：
1：在计算步长时我做错了什么，或者这个例子只是病态的，如果是这样，为什么？
2：如何实际计算提供稳定性的最大步长？

1个回答

使用精确解来估计 $\lambda$ 可能无法给出适用于数值方案的时间步长。当地估计 $\lambda$ 是

λ = 400 t u

$\lambda = 400 t u$ 如果你尝试这一步

h = min (2 / 20, 2 / (400 t u + 10^{- 12}))

$h = \min(2/20, 2/(400 t u+10^{-12}))$ 它是稳定的，但当然准确性很差。这是一个代码

import numpy as np
import matplotlib.pyplot as plt

t = 0.0
u = 1.0
tdata, udata = [], []
tdata.append(t); udata.append(u)
while t<3.0:
    dt = np.min([2.0/20.0, 2.0/(400*t*u+1.0e-12)])
    rhs = - 200.0*t*u**2
    u = u + dt*rhs
    t += dt
    print("dt, t =", dt, t)
    tdata.append(t); udata.append(u)

tdata = np.array(tdata); udata = np.array(udata)
uexact = 1.0/(1.0 + 100*tdata**2)
plt.plot(tdata,udata,tdata,uexact)
plt.legend(('Forward Euler','Exact'))
plt.xlabel('t'); plt.ylabel('u')
plt.show()

第一步，步长为。第一步之后，数值解仍然是和所以你需要在第二步使用步长，小于 . $h=0.1$ $u=1$

λ = 400 * 0.1 * 1 = 40

$\lambda = 400*0.1*1 = 40$

h \leq 2 / 40 = 0.05

$h \le 2/40 = 0.05$

2 / 20 = 0.1

$2/20=0.1$

这个特定问题还有另一个问题。从开始，解决方案应保留在中。如果由于某种原因它变为负数，则正向 euler 会将其带到，就像一样。ODE 本身具有这种行为。所以对于这个问题，如果精确解应该是正的，则应该选择时间步长，使数值解保持为正。我们有 iff 而且这个条件也符合线性稳定性分析的限制。 $u=1$ $[0,1]$ $-\infty$ $h=2/20$

u_{n + 1} = u_{n} - 200 t_{n} u_{n}^{2} h

$u_{n+1} = u_n - 200 t_n u_n^2 h$

u_{n} \geq 0

$u_n \ge 0$

⟹

$\implies$

u_{n + 1} \geq 0

$u_{n+1} \ge 0$

h \leq \frac{1}{200 t_{n} u_{n}}

$h \le \frac{1}{200 t_n u_n}$

其它你可能感兴趣的问题

上一篇在 FEM 模型中定义电流密度 (MATLAB) 下一篇如何确定图中是否存在闭环路径？