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

Question

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

但是，我的代码给出的值与我的书本或 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))

Answer 1

你的方程一般有解

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

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

您可以强制积分器在其内部步骤中更加精确，而无需通过设置误差阈值来减小输出步长 dt

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) ]

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

Why is my code using 4th Runge-Kutta isn't giving me the expected values?

python

numpy

matplotlib

runge-kutta