为什么我使用 4th Runge-Kutta 的代码没有给我预期的值？答案

【问题标题】：Why is my code using 4th Runge-Kutta isn't giving me the expected values?为什么我使用 4th Runge-Kutta 的代码没有给我预期的值？
【发布时间】：2017-04-18 19:56:46
【问题描述】：

我在试图理解我的代码有什么问题时遇到了一些麻烦，任何帮助都会非常有帮助。我想解这个简单的方程

但是，我的代码给出的值与我的书本或 wolfram 的值不匹配，因为 y 随着 x 的增长而上升。

import matplotlib.pyplot as plt
from numpy import exp
from scipy.integrate import ode

# initial values 
y0, t0 = [1.0], 0.0

def f(t, y):
    f = [3.0*y[0] - 4.0/exp(t)]
    return f

# initialize the 4th order Runge-Kutta solver
r = ode(f).set_integrator('dopri5')
r.set_initial_value(y0, t0)

t1 = 10
dt = 0.1
x, y = [], []
while r.successful() and r.t < t1:
    x.append(r.t+dt); y.append(r.integrate(r.t+dt))
    print(r.t+dt, r.integrate(r.t+dt))

【问题讨论】：

试试更小的dt。你会注意到在它爆炸之前你走得更远。当你需要在精确解相当平滑的情况下采取过小的步骤时，你应该怀疑你的 ODE 是“僵硬的”：en.wikipedia.org/wiki/Stiff_equation
哦，我明白了！使用较小的 dt 方程可以正常工作，现在我明白了，我以为我做错了什么，因为这是我第一次使用这个 ode 函数。好的，谢谢。
@klimaat。您能否发表您的评论作为答案，以便 OP 可以接受它并且可以将问题标记为已回答？此外，你们都会从中获得一点声誉，这很好。

标签： python numpy matplotlib runge-kutta

【解决方案1】：

你的方程一般都有解

y(x) = (y0-1)*exp(3*x) + exp(-x)

由于初始条件的选择，精确解不包含第一项的增长分量。然而，由于离散化和浮点误差引起的小扰动将在增长项中产生非零系数。现在，在积分区间结束时，该随机系数乘以 exp(3*10)=1.107e+13，这将放大 1e-7 大小的小离散化误差对大小 1e+6 的结果的贡献，正如运行原始代码时所观察到的那样。

您可以通过设置错误阈值（如

r = ode(f).set_integrator('dopri5', atol=1e-16, rtol=1e-20)

但是，您无法完全避免结果的恶化，因为大小为1e-16 的浮点错误被放大到大小为1e-3 的全局错误贡献中。

此外，您应该注意到，每次调用r.integrate(r.t+dt) 都会将积分器推进dt，以便存储的数组和打印的值处于同步状态。如果您只想打印积分器的当前状态，请使用

    print(r.t,r.y,yexact(r.t,y0))

最后一个是与确切的解决方案进行比较，如前所述，

def yexact(x,y0): 
    return [ (y0[0]-1)*exp(3*x)+exp(-x) ]

【讨论】：