计算科学 - 变系数线性常微分方程IVP的数值解爆破 - 吾爱随笔录

变系数线性常微分方程IVP的数值解爆破

计算科学颂稳定时间积分

2021-12-22 13:10:41

Cross 发布在 Mathematica.SE 中，我将在此处尝试以更一般的方式重新表述它。

我的一个朋友向我展示了这个具有可变系数的线性常微分方程 (ODE) 的初始值问题 (IVP)：

y^{″} (x) = (x^{2} - 1) y (x)

$y''(x)=\left(x^2-1\right) y(x)$

y (0) = 1

$y(0)=1$

y^{'} (0) = 0

$y'(0)=0$

看起来很简单，对吧？实际上它可以解析地解决，解决方案是：

y (x) = e^{- \frac{x^{2}}{2}}

$y(x)=e^{-\frac{x^2}{2}}$

但是当我试图用经典的龙格-库塔方法来解决它时，数值解就崩溃了：

这里使用的步长为 0.001。

为什么会这样？如果我选择了不正确的方法，那么什么是合适的？

3个回答

这是二阶 ODE，因此有两种解决方案。第二种解决方案是不稳定的，因为它要么趋于正无穷，要么趋于负无穷 $x$ 变大。为了 $x$ 接近于零，解显然是稳定的。即使这个解是你上面展示的解析解，在数值积分算法中，由于舍入误差，第二个解的一个分量会潜入。最终，这个分量会增长并主导解。

我相信所有用于求解初始值 ODE 的标准积分算法都会表现出这种行为。

您的分析解决方案假设 $y>=0$ 为了 $x>0$ . 然而，对于初始条件 $y<0$ , $y(x) = -e^{-x^2/2}$ . 详细分析该系统的稳定性可能是题外话，并且更适合数学站点。

解渐近收敛到 $0$ 成倍增长。在数值上，大多数集成方法都会遇到困难。小值周围的数值噪声会导致求解器采取过大的步长，从而导致 $y(x)$ 去消极。解决此问题的一种方法是使用能够智能地将解决方案组件约束为非负的求解器。正确地做到这一点并非易事——参见Shampine 等人。2005 年。Matlab 的 ODE 求解器，例如，通过可以使用设置的输出属性ode45具有此功能。我不认为数学有这样的选择。这个Mathematica.StackExchange 帖子提供了几种可以尝试添加它的方法。'NonNegative' odesetNDSolve

这是 Matlab 中的一个简单示例：

f = @(x,y)[y(2);(x^2-1)*y(1)];
opts = odeset('NonNegative',1); % Keep y(1) >= 0
[x,y] = ode45(f,[0 10],[1 0],opts);

绘图，你会得到这样的东西：

模拟输出图

完整的解决方案是

y (x) = C \cdot e^{- x^{2} / 2} + D \cdot e^{- x^{2} / 2} \cdot \int_{0}^{x} e^{s^{2}} d s

$y(x)=C·e^{-x^2/2}+D·e^{-x^2/2}·\int_0^x e^{s^2}ds$ 第一个是偶函数并且有界，第二个是奇函数，具有不可进一步约化的积分并且像爆炸一样

\sim x \cdot e^{x^{2} / 2}

$\sim x·e^{x^2/2}$ .

其它你可能感兴趣的问题

上一篇使用可变长度的布尔变量序列作为索引来划分数组下一篇Galerkin FEM和Discontinuous Galerkin中保守离散化的意义