我对精度为 25 到 30 位有效数字的常微分方程 (ODE) 的数值解感兴趣（我的意思是在每个点找到实际解的 25 到 30 位有效数字）。
我的第一个猜测是多精度或任意精度算术。
但是根据这个post，有多个精度包可以在没有多个精度的情况下进行计算（？！）。这可能并不像听起来那么矛盾。然而，它带来了一个严重的问题。
那么是否存在任何可以进行可靠的多精度计算的多精度包（在 Fortran、C、C++、Python、Java 或其他中）？
【问题讨论】：
标签： precision differential-equations
出于演示目的，让我们保持谦虚，并在圆方程 u'=v; v'=-u 的实现中瞄准 20 个正确的数字，这使得 u=sin(t) 对应于 (u0,v0)=(0,1)。让我们在区间 [0,1] 上整合它。
这有一个 L=1 的 Lipschitz 常数，这使得争论步长更容易。使用经典的 Runge-(Heun)-Kutta 4 阶方法，全局误差预计为 h^4 的一小部分，因此要在点后获得大约 20 位数字，我们需要 h=1e-5，这也意味着 100 000 步骤在区间内。每个数字的最后一位都有一些浮点错误，因此在整个区间内，可能有 5 位数字是不可靠的。因此工作精度应为25位或更多
from mpmath import mp
mp.dps = 25
def RK4(t0, y0, tf, N):
h = (tf-t0)/N
u,v = y0;
for n in range(N):
k1u, k1v = v, -u;
k2u, k2v = v+0.5*h*k1v, -(u+0.5*h*k1u)
k3u, k3v = v+0.5*h*k2v, -(u+0.5*h*k2u)
k4u, k4v = v+h*k3v, -(u+h*k3u)
u, v = u+h*(k1u+2*k2u+2*k3u+k4u)/6, v+h*(k1v+2*k2v+2*k3v+k4v)/6
return u,v
u,v = mp.mpf(0), mp.mpf(1)
t, dt = mp.mpf(0), mp.mpf(1)/10
for k in range(10):
u1,v1 = RK4(t,[u,v],t+dt,10**3)
u,v = RK4(t,[u,v],t+dt,10**4)
t = t+dt
# error estimate via extrapolation
# u = uexact + uerr, u1 = uexact + 10**4*uerr
uerr, verr = (u1-u)/9999, (v1-v)/9999
with mp.extradps(5): si, co = mp.sin(t), mp.cos(t)
print("t=%s"%t)
print("u=%30s, sin(t)=%30s, step error=%30s"%(u, si, uerr))
print("v=%30s, cos(t)=%30s, step error=%30s"%(v, co, verr))
给出如下结果列表
t=0.1
u= 0.09983341664682815230680593, sin(t)= 0.0998334166468281523068142, step error=-8.291772901773692681703852e-24
v= 0.9950041652780257660955634, cos(t)= 0.995004165278025766095562, step error=8.312005953404312840899088e-25
t=0.2
u= 0.1986693307950612154593963, sin(t)= 0.1986693307950612154594126, step error=-8.167351203006133449374218e-24
v= 0.9800665778412416311242006, cos(t)= 0.9800665778412416311241965, step error=1.654857583109387499562722e-24
t=0.3
u= 0.2955202066613395751052971, sin(t)= 0.2955202066613395751053207, step error= -7.9613619853238077687103e-24
v= 0.9553364891256060196423186, cos(t)= 0.9553364891256060196423102, step error= 2.46199127828304457142136e-24
t=0.4
u= 0.3894183423086504916662812, sin(t)= 0.3894183423086504916663118, step error=-7.675762243562770001241619e-24
v= 0.9210609940028850827985405, cos(t)= 0.9210609940028850827985267, step error=3.244534570708288054285545e-24
t=0.5
u= 0.479425538604203000273252, sin(t)= 0.4794255386042030002732879, step error=-7.313560501818737799847937e-24
v= 0.8775825618903727161163026, cos(t)= 0.8775825618903727161162816, step error=3.994583217701846864434659e-24
t=0.6
u= 0.5646424733950353572009052, sin(t)= 0.5646424733950353572009454, step error=-6.878230619815552707865177e-24
v= 0.8253356149096782972409814, cos(t)= 0.8253356149096782972409525, step error=4.704815938714571792025104e-24
t=0.7
u= 0.6442176872376910536725712, sin(t)= 0.6442176872376910536726143, step error=-6.374164848843117446876767e-24
v= 0.7648421872844884262558976, cos(t)= 0.76484218728448842625586, step error=5.368009690718806448531807e-24
t=0.8
u= 0.7173560908995227616271302, sin(t)= 0.7173560908995227616271746, step error=-5.806451504735070057940143e-24
v= 0.6967067093471654209207975, cos(t)= 0.69670670934716542092075, step error=5.977494663044775070749797e-24
t=0.9
u= 0.7833269096274833884613378, sin(t)= 0.7833269096274833884613823, step error=-5.180615155476414729415392e-24
v= 0.621609968270664456484775, cos(t)= 0.6216099682706644564847162, step error=6.527225370323768115169452e-24
t=1.0
u= 0.8414709848078965066524596, sin(t)= 0.8414709848078965066525023, step error=-4.503113625506337452188567e-24
v= 0.540302305868139717401006, cos(t)= 0.5403023058681397174009366, step error=7.011939642161084587215687e-24
可以看到错误仅限于最后 4 位，正如预期的那样。每个子区间的累积方法误差约为1e-23，因此总共小于1e-22，与预期的浮点噪声大致相同，它在准随机舍入或算术运算的总和中会影响最多 5 位数，平均约为一半。
高阶方法将允许更大的步长，从而减少步长，这不仅减少了计算时间，而且还减少了舍入误差的累积。因此，工作精度需要更少的额外数字，这也（略微）减少了计算时间。
对于像这样的繁重计算，应该使用编译语言。强类型化和从垃圾收集到主动内存管理的转换将使代码看起来不那么容易阅读，就像在How are the trigonometric functions tested in the GNU C Library?中引用的代码一样
【讨论】：