我对普通微分方程（ODE）的数值解感兴趣，其精度为25到30个有效数字（我的意思是在每个点上找到实际解的25到30个有效数字）。
我的第一个猜测是多精度或任意精度算术。
然而，根据此post，存在多个精度包，它们执行计算时没有多个精度（？！）。这可能听起来并不矛盾。但是，这带来了严重的问题。
那么，是否存在任何可以执行可靠的多精度计算的多精度软件包（在Fortran，C，C ++，Python，Java或任何其他程序中？
u'=v; v'=-u的u=sin(t)的圆方程(u0,v0)=(0,1)的实现中，争取20个正确的数字。让我们在间隔[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?]中引用的代码一样