有数千篇论文和数百个代码使用五阶或更高阶的 Runge-Kutta 方法。请注意，MATLAB 中最常用的显式积分器是 ODE45，它使用 5 阶 Runge-Kutta 方法推进解决方案。

广泛使用的高阶龙格-库塔方法示例

根据Google Scholar ，Dormand & Prince 给出 5 阶方法的论文有超过 1700 次引用。其中大多数是使用他们的方法解决某些问题的论文。 Cash-Karp 方法论文有超过 400 次引用。也许最广泛使用的高于 5 的阶数方法是 Prince-Dormand 的 8 阶方法，它在 Google Scholar 上有超过 400 次引用。我可以举很多其他的例子；请记住，许多（如果不是大多数）使用这些方法的人从不引用论文。

另请注意，高阶外推法和延迟校正方法是龙格-库塔法。

高阶方法和舍入误差

如果您的准确性受到舍入误差的限制，那么您应该使用高阶方法。这是因为高阶方法需要更少的步骤（和更少的函数评估，即使每个步骤有更多的评估），所以它们提交的舍入错误更少。您可以通过简单的实验轻松验证这一点；对于数值分析的第一门课程来说，这是一个很好的家庭作业问题。

十阶方法在双精度算术中非常有用。相反，如果我们只有欧拉方法，那么舍入误差将是一个主要问题，我们需要非常高精度的浮点数来解决高阶求解器可以解决的许多问题。

高阶方法可以同样稳定

@RichardZhang 引用了第二个 Dahlquist 障碍，但这仅适用于多步方法。这里发布的问题是关于 Runge-Kutta 方法的，每个订单都有 Runge-Kutta 方法，不仅$A$- 稳定，但也$B$-stable（对一些非线性问题有用的稳定性属性）。要了解这些方法，请参阅例如 Hairer & Wanner 的文本。

天体力学中的高阶方法

你问

当应用于在极端时间尺度上具有高度振荡解的方程时，这种高阶方法通常不是首选吗？

你完全正确！天体力学就是一个典型的例子。我不是那个领域的专家。但是，例如，这篇论文比较了天体力学的方法，甚至没有考虑低于 5 的阶数。它得出的结论是 11 阶或 12 阶的方法通常是最有效的（8 阶的 Prince-Dormand 方法也经常非常高效的）。

基准设置

在 Julia 软件DifferentialEquations.jl中，我们实现了许多高阶方法，包括 Feagin 方法。您可以在我们的方法列表中看到它，然后还有很多其他方法可以用作提供的画面。因为所有这些方法都放在一起，您可以轻松地在它们之间进行基准测试。您可以在此处查看我在网上的基准测试，并看到对许多不同算法进行基准测试非常简单。因此，如果您想花几分钟时间运行基准测试，那就去做吧。以下是结果的摘要。

首先需要注意的是，如果您查看每个基准测试，您会发现我们的DP5（Dormand-Prince Order 5）和DP8方法比 Hairer Fortran 代码（dopri5和dop853）更快，因此这些实现得到了很好的优化. 这些表明，正如另一个线程中所指出的，过度使用 Dormand-Prince 方法是因为这些方法已经编写好了，而不是因为它们仍然是最好的。因此，最优化实现之间的真正比较是在 Tsitorous 方法、Verner 方法和来自 DifferentialEquations.jl 的 Feagin 方法之间进行比较。

结果

一般来说，阶数高于 7 的方法具有额外的计算成本，考虑到所选择的容差，阶数通常不会超过该成本。原因之一是低阶方法的系数选择更加优化（它们具有小的“原则截断误差系数”，当您不是渐近地小时，这更重要）。您可以看到，在像这里这样的许多问题中，Verner Efficient 7 和 8 方法都做得非常好，但是像 Verner Efficient 9 这样的方法可以有一个较低的斜率。这是因为高阶的“增益”在较低的容差下复合，所以总是存在高阶方法更有效的容差。

然而，问题是，有多低？在一个优化良好的实现中，由于两个原因，它变得相当低。第一个原因是因为低阶方法实现了称为 FSAL 的东西（第一个与最后一个相同）。此属性意味着低阶方法在下一步中重复使用上一步中的函数评估，因此有效地减少了一次函数评估。如果使用得当，那么像 5 阶方法（Tsitorous 或 Dormand-Prince）之类的方法实际上会进行 5 次函数评估，而不是表格建议的 6 次。Verner 6 方法也是如此。

另一个原因是由于插值。使用非常高阶方法的一个原因是减少步骤并简单地插入中间值。但是，为了获得中间值，插值函数可能需要比用于采取步骤更多的函数评估。如果您查看 Verner 方法，对于 8 阶方法，需要 8 次额外的函数评估才能获得 8 阶插值。很多时候，低阶方法提供“免费”插值，例如大多数 5 阶方法具有免费的 4 阶插值（没有额外的函数评估）。因此，这意味着如果您需要中间值（如果您使用高阶方法，则需要中间值来获得良好的绘图），会有一些额外的隐藏成本。考虑到这些插值对于事件处理和求解延迟微分方程非常重要这一事实，您就会明白为什么额外的插值成本会考虑在内。

那么Feagin方法呢？

因此，您会看到基准测试中可疑地缺少 Feagin 方法。它们很好，收敛测试适用于任意精度数等，但要真正让它们做得好，你需要要求一些非常荒谬的低容差。例如，我在未发表的基准测试中发现Feagin14，Vern9在1e-30. 对于混沌动力学的应用（如昴宿星或三体天体物理学问题），由于敏感的依赖性（混沌系统中的错误快速复合），您可能需要这种准确度。但是，大多数人可能都在使用双精度浮点数进行计算，而且我还没有找到一个基准，它们在这个容差范围内表现出色。

这是这个 100 维 ODE的工作精度效率图，展示了这一点：

这是 Pleiades ODE的工作精度效率图：

有了这些，我们在测试具有 Float64 范围内容差适应性的 Feagin 方法时没有找到好的结果。问题可能是误差估计的稳定性。请注意，这是通过 PI 自适应来完成的，以稳定每个方法的步进。

此外，没有与 Feagin 方法一起使用的插值法（Edit 2020: See Feagin's remark that one now exists, but is not available online。但请注意，这对基准测试没有影响，因为基准测试不使用插值法）。所以我所做的只是简单地对它们进行三阶 Hermite 插值，这样就可以存在（而且效果非常好）。但是，如果没有标准的插值函数，您可以使用递归 Hermite 方法（使用此插值得到中点，然后进行 5 阶插值等）来获得高阶插值，但这非常昂贵且结果插值不一定有一个低原则截断误差项（所以只有在dt真的很小，这与我们想要的情况完全相反！）。因此，如果您需要一个非常好的插值来匹配您的准确性，您至少需要回到类似Vern9.

关于外推的注意事项

请注意，外推方法只是生成任意阶龙格-库塔方法的简单算法。然而，对于它们的顺序，它们采取的步骤比必要的要多，并且具有较高的原则截断误差系数，因此它们不如给定顺序的良好优化的 RK 方法有效。但是鉴于前面的分析，这意味着存在一个容差极低的域，这些方法将比“已知”的 RK 方法做得更好。但是在我跑过的每一个基准测试中，似乎我都没有得到那么低。

关于稳定性的说明

选择确实与稳定性问题无关。事实上，如果您查看 DifferentialEquations.jl 表格（您可以只plot(tab)针对稳定区域），您会发现大多数方法都有可疑的相似稳定区域。这其实是一种选择。通常在推导方法时，作者通常会做以下事情：

1. 找到最低的原则截断误差系数（即下阶项的系数）
2. 受订单约束
3. 并使稳定区域接近 Dormand-Prince Order 5 方法的稳定区域。

为什么是最后一个条件？好吧，因为这种方法在 PI 控制的自适应步长选择完成的方式下往往总是稳定的，所以它是“足够好”稳定区域的一个很好的标准。因此，稳定区域都趋于相似并非巧合。

结论

每种方法的选择都需要权衡取舍。最高阶的 RK 方法在较低的容差下效率不高，这既是因为更难优化系数的选择，也是因为函数评估的数量复合（并且在涉及插值时增长更快）。但是，如果容差变得足够低，它们就会胜出，但所需的容差可能远低于“标准”应用程序（即实际上只适用于混沌系统）。

附录：基准代码

using OrdinaryDiffEq, DiffEqDevTools, Plots
using Random
Random.seed!(123)
gr()
# 2D Linear ODE
function f(du,u,p,t)
  @inbounds for i in eachindex(u)
    du[i] = 1.01*u[i]
  end
end
function f_analytic(u₀,p,t)
  u₀*exp(1.01*t)
end
tspan = (0.0,10.0)
prob = ODEProblem(ODEFunction(f,analytic=f_analytic),rand(100,100),tspan)

abstols = 1.0 ./ 10.0 .^ (3:13)
reltols = 1.0 ./ 10.0 .^ (0:10);

setups = [Dict(:alg=>Tsit5())
          Dict(:alg=>Vern7())
          Dict(:alg=>Vern8())
          Dict(:alg=>Feagin10())
          Dict(:alg=>Feagin12())
          Dict(:alg=>Feagin14())
          Dict(:alg=>Vern9())]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;save_everystep=false,numruns=100)
plot(wp)
savefig("feagin1.png")

f = (du,u,p,t) -> begin
  @inbounds begin
  x = view(u,1:7)   # x
  y = view(u,8:14)  # y
  v = view(u,15:21) # x′
  w = view(u,22:28) # y′
  du[1:7] .= v
  du[8:14].= w
  for i in 14:28
    du[i] = zero(u[1])
  end
  for i=1:7,j=1:7
    if i != j
      r = ((x[i]-x[j])^2 + (y[i] - y[j])^2)^(3/2)
      du[14+i] += j*(x[j] - x[i])/r
      du[21+i] += j*(y[j] - y[i])/r
    end
  end
  end
end

prob = ODEProblem(f,[3.0,3.0,-1.0,-3.0,2.0,-2.0,2.0,3.0,-3.0,2.0,0,0,-4.0,4.0,0,0,0,0,0,1.75,-1.5,0,0,0,-1.25,1,0,0],(0.0,3.0))

abstols = 1.0 ./ 10.0 .^ (9:12)
reltols = 1.0 ./ 10.0 .^ (6:9);

sol = solve(prob,Vern8(),abstol=1/10^12,reltol=1/10^10,maxiters=1000000)
test_sol = TestSolution(sol);

setups = [Dict(:alg=>Tsit5())
          Dict(:alg=>Vern7())
          Dict(:alg=>Vern8())
          Dict(:alg=>Feagin10())
          Dict(:alg=>Feagin12())
          Dict(:alg=>Feagin14())
          Dict(:alg=>Vern9())]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,save_everystep=false,numruns=100,maxiters=1000)
plot(wp)
savefig("feagin2.png")

只要您使用标准的双精度浮点运算，就不需要非常高阶的方法来以合理的步数获得高精度的解决方案。在实践中，我发现解决方案的准确性通常被双精度浮点表示限制为 1.0e-16 的相对误差，而不是使用 RKF45 采取的步骤的数量/长度。

如果您切换到一些高于双精度的浮点算术方案，那么使用 10 阶方法可能是值得的。

好吧，在回答有关使用高阶 RK 方法的问题时，我有几件事要补充到讨论中：

1. 例如，绘制“效率图”对于应用于给定初始值问题的几种不同方法（使用各种步长）通常很有帮助。在这样的图表上，可以绘制（对于每种方法）精度（通常根据获得的有效数字的数量来衡量）如何随log10用于在特定间隔内积分的数字函数/导数调用而变化。

   在这样的图上，可以同时绘制几种不同的方法，以便很好地了解哪种方法性能最好（或为给定的所需精度提供最有效的集成）。在这样的图上，每种方法的直线斜率大约等于方法的阶数（对于那些以截断误差为主的步长）。

   当然，最终，随着步长变得非常小，函数/导数调用的数量变得非常大，舍入误差占主导地位，曲线或多或少地向下弯曲。在效率图 上，对于所采取的相当大的步长，舍入误差是不可见的。

对于该图，要解决的问题是偏心率为 0.4 的二体问题，并且要比较的方法是一些 4、6、8、10、12 阶的 Explicit RK 方法。可以观察到，如果你只想要 4 -6 位数的准确度，那么您可能希望使用哪种方法并没有太大的区别。如果您想要更高的准确性，那么您采用哪种方法变得更加重要。

图中显示的方法包括 RK4、Rk6B（Butcher 的 6 阶方法）、RK8CV（Cooper 和 Verner 的 8 阶方法）、RK10H（Hairer 的 10 阶方法）和 RK12（我的 12 阶方法，它恰好有一个嵌入的 10 阶方法，以便您可以估计局部截断误差）。没有显示我的 14 阶方法，除非您需要极高的精度，否则它不是特别有效），正如上面其他人所指出的。

2. 众所周知，我的 RK12(10) 确实有能力提供 7 阶的插值结果。您仍然可以在每个步骤结束时获得 12 阶结果，此外，如果您需要所谓的“密集输出” "，您可以在步骤中的任何时候使用 7 阶结果。描述该方法的最佳方式是让我向您发送一份代码副本。我宁愿不在网站上发布代码。这样，我可以确保任何更新等都可以分发给感兴趣的人。

   只需将任何请求发送至 feagin@uhcl.edu（我在休斯顿大学 - 明湖教授计算机科学）。