python - 需要更好地了解 rtol、atol 在 scipy.integrate.odeint 中的工作方式

Question

这里 scipy.integrate.odeint 被调用了六个不同的标准 ode 问题，rtol = atol 来自 1E-06 到 1E-13。我查看了所有较大公差减去最小公差的结果之间的最大差异，以获得某种“错误”表示。我很好奇为什么对于给定的公差，一个问题 (D5) 给出的错误比另一个问题 (C1) 严重一百万倍，即使步骤数的范围相当小(在 10 的因数内)。

脚本中给出了颂歌问题的引用。所有问题都相当规范化，所以我以类似的方式处理 rtol 和 atol。

重申一下 - 我的问题是，为什么不同问题之间的误差几乎相差 1E+06，尽管误差随公差而变化。当然，C1 是“最软的”，D5 在“近日点”处有戏剧性的峰值，但我认为该例程会在内部调整步长，以便误差相似。

编辑:我添加了“错误”的时间演化，这可能会说明一些问题。

# FROM: "Comparing Numerical Methods for Ordinary Differential Equations"
# T.E. Hull, W.H. Enright, B.M. Fellen and A.E. Sedgwidh
# SIAM J. Numer. Anal. vol 9, no 4, December 1972, pp: 603-637

def deriv_B1(y, x):
    return [2.*(y[0]-y[0]*y[1]), -(y[1]-y[0]*y[1])] # "growth of two conflicting populations"

def deriv_B4(y, x):
    A = 1./np.sqrt(y[0]**2 + y[1]**2)
    return [-y[1] - A*y[0]*y[2],  y[0] - A*y[1]*y[2],  A*y[0]]  # "integral surface of a torus"

def deriv_C1(y, x):
    return [-y[0]] + [y[i]-y[i+1] for i in range(8)] + [y[8]] # a radioactive decay chain

def deriv_D1toD5(y, x):
    A = -(y[0]**2 + y[1]**2)**-1.5
    return [y[2],  y[3],  A*y[0],  A*y[1]] # dimensionless orbit equation

deriv_D1, deriv_D5 = deriv_D1toD5, deriv_D1toD5

def deriv_E1(y, x):
    return [y[1], -(y[1]/(x+1.0) + (1.0 - 0.25/(x+1.0)**2)*y[0])] # derived from Bessel's equation of order 1/2

def deriv_E3(y, x):
    return [y[1], y[0]**3/6.0 - y[0] + 2.0*np.sin(2.78535*x)] # derived from Duffing's equation

import numpy as np
from scipy.integrate import odeint as ODEint
import matplotlib.pyplot as plt
import timeit

y0_B1 = [1.0, 3.0]
y0_B4 = [3.0, 0.0, 0.0]
y0_C1 = [1.0] + [0.0 for i in range(9)]
ep1, ep5 = 0.1, 0.9
y0_D1 = [1.0-ep1, 0.0, 0.0, np.sqrt((1.0+ep1)/(1.0-ep1))]
y0_D5 = [1.0-ep5, 0.0, 0.0, np.sqrt((1.0+ep5)/(1.0-ep5))]
y0_E1 = [0.6713968071418030, 0.09540051444747446] # J(1/2, 1), Jprime(1/2, 1)
y0_E3 = [0.0, 0.0]

x  = np.linspace(0, 20, 51)
xa = np.linspace(0, 20, 2001)

derivs = [deriv_B1, deriv_B4, deriv_C1, deriv_D1, deriv_D5, deriv_E3]
names  = ["deriv_B1", "deriv_B4", "deriv_C1", "deriv_D1", "deriv_D5", "deriv_E3"]
y0s    = [y0_B1, y0_B4, y0_C1, y0_D1, y0_D5, y0_E3]

timeit_dict, answer_dict, info_dict = dict(), dict(), dict()

ntimes = 10
tols   = [10.**-i for i in range(6, 14)]

def F():           # low density of time points, no output for speed test
    ODEint(deriv, y0, x, rtol=tol, atol=tol)
def Fa():           # hight density of time points, full output for plotting
    return ODEint(deriv, y0, xa, rtol=tol, atol=tol, full_output=True)

for deriv, y0, name in zip(derivs, y0s, names):
    timez = [timeit.timeit(F, number=ntimes)/float(ntimes) for tol in tols]
    timeit_dict[name] = timez
    alist, dlist = zip(*[Fa() for tol in tols])
    answer_dict[name] = np.array([a.T for a in alist])
    info_dict[name] = dlist

plt.figure(figsize=[10,6])

for i, name in enumerate(names):
    plt.subplot(2, 3, i+1)
    for thing in answer_dict[name][-1]:
        plt.plot(xa, thing)
    plt.title(name[-2:], fontsize=16)
plt.show()

plt.figure(figsize=[10, 8])
for i, name in enumerate(names):
    plt.subplot(2,3,i+1)
    a = answer_dict[name]
    a13, a10, a8 = a[-1], a[-4], a[-6]
    d10 = np.abs(a10-a13).max(axis=0)
    d8  = np.abs(a8 -a13).max(axis=0)
    plt.plot(xa, d10, label="tol(1E-10)-tol(1E-13)")
    plt.plot(xa, d8,  label="tol(1E-08)-tol(1E-13)")
    plt.yscale('log')
    plt.ylim(1E-11, 1E-03)
    plt.title(name[-2:], fontsize=16)
    if i==3:
        plt.text(3, 1E-10, "1E-10 - 1E-13", fontsize=14)
        plt.text(2, 2E-05, "1E-08 - 1E-13", fontsize=14)
plt.show()

fs = 16
plt.figure(figsize=[12,6])

plt.subplot(1,3,1)
for name in names:
    plt.plot(tols, timeit_dict[name])
plt.title("timing results", fontsize=16)
plt.xscale('log')
plt.yscale('log')
plt.text(1E-09, 5E-02, "D5", fontsize=fs)
plt.text(1E-09, 4.5E-03, "C1", fontsize=fs)

plt.subplot(1,3,2)
for name in names:
    a = answer_dict[name]
    e = a[:-1] - a[-1]
    em = [np.abs(thing).max() for thing in e]
    plt.plot(tols[:-1], em)
plt.title("max difference from smallest tol", fontsize=16)
plt.xscale('log')
plt.yscale('log')
plt.xlim(min(tols), max(tols))
plt.text(1E-09, 3E-03, "D5", fontsize=fs)
plt.text(1E-09, 8E-11, "C1", fontsize=fs)

plt.subplot(1,3,3)
for name in names:
    nsteps = [d['nst'][-1] for d in info_dict[name]]
    plt.plot(tols, nsteps, label=name[-2:])
plt.title("number of steps", fontsize=16)
plt.xscale('log')
plt.yscale('log')
plt.ylim(3E+01, 3E+03)
plt.legend(loc="upper right", shadow=False, fontsize="large")
plt.text(2E-12, 2.3E+03, "D5", fontsize=fs)
plt.text(2E-12, 1.5E+02, "C1", fontsize=fs)

plt.show()

Answer 1

自从我发布问题后，我学到了更多。不能只把每一步的数值精度乘以步数，就希望得到整体的精度。

如果解决方案出现分歧(附近的起点导致路径随着时间的推移变得越来越远)，那么数值误差就会被放大。每个问题都会有所不同 - 一切都应该如此。

赫尔等人。是学习 ODE 求解器的一个很好的起点。 (问题中显示的问题的来源)

"Comparing Numerical Methods for Ordinary Differential Equations" T.E. Hull, W.H. Enright, B.M. Fellen and A.E. Sedgwidh SIAM J. Numer. Anal. vol 9, no 4, December 1972, pp: 603-637