需要更好地了解rtol，atol在scipy.integrate.odeint中的工作方式 [英] need to understand better how rtol, atol work in scipy.integrate.odeint
问题描述
在scipy.integrate.odeint被调用时有六个不同的标准ode问题，其中rtol = atol从1E-06到1E-13.我已经查看了所有较大容差的结果之间的最大差值减去最小容差的结果之间的最大差值，以得到某种形式的错误"表示.我很好奇为什么在给定的公差下，即使步数的范围非常狭窄(在10倍以内)，一个问题(D5)的错误也比另一个问题(C1)的错误要严重一百万倍.
Here scipy.integrate.odeint is called with six different standard ode problems with rtol = atol from 1E-06 to 1E-13. I've looked at the max difference between the results at all larger tolerances minus those of the smallest, to get some kind of representation of "error". I'm curious why, for a given tolerance, one problem (D5) gives errors a million times worse than another problem (C1), even though the range in number of steps is fairly tight (within a factor of 10).
在脚本中给出了对ode问题的引用.所有问题都已很好地标准化，因此我将类似地对待rtol和atol.
The citation for the ode problems is given in the script. All problems are fairly well normalized so I'm treating rtol and atol similarly.
重申一下-我的问题是，为什么误差在不同问题之间的变化幅度几乎为1E+06，尽管误差会随着公差的变化而变化.当然C1是最软"的，而D5在"perihelion"处具有惊人的峰值，但我当时认为例程会在内部调整步长，以使误差相似.
To reiterate - my question is why the errors vary by a factor of almost 1E+06 between different problems, though the errors scale with tolerance. Of course C1 is the "softest" and D5 has the dramatic peaks at "perihelion" but I was thinking that the routine would adjust the step sizes internally so that the errors would be similar.
编辑:我添加了错误"的时间演变，这可能会引起一些启发.
I've added the time evolution of the "errors" which may shed some light.
# 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()
推荐答案
自从我发布问题以来，我已经学到了更多.不能只将每步的数值精度乘以步数，并希望获得整体精度.
Since I posted the question, I've learned more. One can't just multiply the numerical accuracy per step by the number of steps, and hope to get the overall accuracy.
如果解决方案有所不同(附近的起点导致路径随着时间的推移变得越来越远)，则数字误差可能会放大.每个问题都会有所不同-一切都应如此.
If solutions diverge (nearby starting points lead to paths which become much farther apart over time) then numerical errors can become amplified. Every problem will be different - all is as it should be.
Hull等.是学习ODE求解器的理想起点. (问题中显示的问题的来源)
Hull et al. is a great place to start when learning about ODE solvers. (the source for the problems shown in the question)
比较常微分方程的数值方法" T.E.赫尔(W.H.) Enright，B.M.费伦和塞奇维德(A.E. Sedgwidh) 暹罗·努默(SIAM J.肛门第9卷，第4期，1972年12月，第603-637页
"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
这篇关于需要更好地了解rtol，atol在scipy.integrate.odeint中的工作方式的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
