今天，我要欢迎Josh Meyer作为本周的客座博主。Josh在MathWorks的文档团队工作，在那里他撰写和维护一些MATLAB数学文档。在这篇文章中，Josh提供了一些关于如何选择使用哪种ODE求解器的建议。交给你了，Josh…

内容

初值问题

MATLAB中有7个常微分方程初值问题求解器:

(注意,ode15i被排除在讨论之外，因为它解决了自己的初值问题:完全隐式的ode形式为$f(t,y,y') = 0$)

要在求解器之间进行选择，首先需要了解为什么对于给定的问题，一个求解器可能比另一个更好。

MATLAB中的ODE求解器都处理初值问题，

$$y' = f \left(t,y \right)$$

其中$y' = dy/dt$。还有一种更普遍的形式，

$$ M(t,y) y' = f \left(t,y \right)$$

其中$M(t,y)$被称为质量矩阵．

从初始条件$y_0$和得到答案的一段时间$(t_0,t_f)$开始，根据求解器的算法，使用前面步骤的结果迭代地得到解。在第一个这样的步骤中，初始条件提供了允许集成继续进行的必要信息。最后的结果是ODE求解器返回一个时间步长为$t_0,t_1，…，t_f$以及每个时间步$y_0,y_1，…，y_f$对应的解。

从理论上讲，这种数值求解技术是可能的，因为微积分基本定理提供了微分方程和积分之间的联系:

$ $ y (t + h) = y (t) + \ int_t ^ {t + h} f \离开(年代,y (s) \右)ds $ $

计算y(t+h)的问题变成了如何近似右边的积分。这就是不同求解器的用武之地。每个不同的求解器都使用不同的数值技术来评估积分，每个求解器都在效率和精度之间进行权衡。

例子:欧拉方法

欧拉方法是一个简单的ODE求解器，但它提供了一个ODE求解器算法中效率和精度之间权衡的例子。假设你想解

$$y' = f(t,y) = 2t$$

在时间跨度$[0,3]$上使用初始条件$y_0 = 0$。欧拉方法的每一步都是用

数组$ $ \开始{}{cl} y_ {n + 1} & =最大+ h f \左(右t_n,推出\)\ \ t_ {n + 1}识别& = t_n + h \{数组}$ $

使用$h=1$，解决方案只需要三个步骤:

数组$ $ \开始{}{cl} y_1 & = y_0 + f \离开(t_0, y_0 \右)= 0 \ \ y_2 & = y_1 + f \离开(t_1, y_1 \右)= 2 \ \ y_3 & = y_2 + f \离开(t_2, y_2 \右)= 6 \{数组}$ $

．.．但它准确吗?

不是真的。这个方程的精确解是

y(t) = t^2

减少步长$h$可以提高一点答案的准确性，但它也需要更多的步骤来实现解决方案。为了看到这一点，下面的代码使用欧拉方法解决了这个问题，并将答案与$h$的几个不同值的解析解进行了比较。

清晰，CLC h = 1;Tspan = [0 3];F = @(t,y) 2*t;Dydt (1) = 0;T (1) = 0;Y = @(t) t.^2;X = linspace(0,tspan(end));图(x,y(x))“t”), ylabel (“y (t)”)举行在而H > 0.1为k = 2: tspan(结束)/ h + 1 dydt (k) = dydt (k - 1) + h * f (t (k - 1), dydt (k - 1));T (k) = T (k-1) + h;结束情节(t, dydt“o”) h = h/2;结束传奇(精确解的，“h = 1”，“h = 0.5”，“h = 0.25”，“h = 0.125”，．.．“位置”，“西北”)标题(多步长欧拉法求解y”=2t）

欧拉方法的改进

使用越来越小的步长并不是一个好主意，因为算法会降低效率。对于任何合理的问题，这样的求解器都将非常缓慢。此外，欧拉方法也有一些固有的问题。由于$y$的斜率只在每个区间的开始处计算一次，因此这个求解器只对常数函数产生精确的答案。也没有办法估计误差，所以求解器需要使用固定的步长。

因此，改进欧拉方法的一种方法是在每一步中更频繁地计算y'。这提供了中间斜率，可以更好地了解函数在每个区间内的情况，使求解器能够为更高阶的问题产生精确的答案。例如，如果您在欧拉方法中添加了对每个区间中间斜率的计算，则结果称为中点规则，可对线性函数进行精确积分:

数组$ $ \开始{}{cl} s_1 & = f (t_n推出)\ \ s_2 & = f \离开(t_n + \压裂{h}{2},推出+ \压裂{h} {2} s_1 \) \ \ y_ {n + 1} & =最大+ hs_2 \ \ t_ {n + 1}识别& = t_n + h \{数组}$ $

如果在每个区间求四次斜率，就会得到经典龙格-库塔算法(或称。RK4)，这是一块数值算法。该算法对三次函数产生精确的积分(如果$f$只是$t$的函数，则$s_2=s_3$，这与辛普森法则交):

数组$ $ \开始{}{cl} s_1 & = f (t_n推出)\ \ s_2 & = f \离开(t_n + \压裂{h}{2},推出+ \压裂{h} {2} s_1 \) \ \ s_3 & = f \离开(t_n + \压裂{h}{2},推出+ \压裂{h} {2} s_2 \) \ \ s_4 & = f \离开(t_n + h,推出+ hs_3 \) \ \ y_ {n + 1} & =最大+ \压裂{h}{6} \离开(s_1 + 2 s_2 + 2 s_3 + s_4 \) \ \ t_ {n + 1}识别& = t_n + h \{数组}$ $

龙格-库塔算法单步执行求解器，因为每一步都只依赖于前一步的结果。数值，ode23，ode23s，ode23t,ode23tb所有算法都采用单步算法。多步算法，比如ode113而且ode15s，使用过去几个步骤的结果。

复杂的ODE求解器，比如MATLAB中的那些，也会估计每一步中的误差，以确定下一步的大小应该有多大。这是对上面使用的固定步长的另一个改进，因为求解器每一步做更多的工作，可以通过采取不同大小的步长进行补偿。用于确定步长的误差估计通常是通过比较两种不同方法的结果来获得的。MATLAB的ODE求解器遵循一个命名约定，该约定揭示了它们使用的方法的信息。数值比较了四阶龙格-库塔方法和五阶龙格-库塔方法的结果，确定误差。同样的,ode23使用了二阶和三阶龙格-库塔比较。一般来说，数值越小odeNN，解算器的容错越宽松。

那么，这就不足为奇了数值得到了我们之前用欧拉方法解出的方程的一个非常精确的答案。数值是MATLAB的通用ODE求解器，它是你应该用于大多数问题的第一个求解器。

Y = @(t) t.^2;X = linspace(0,3);图(x,y(x))“t”), ylabel (“y (t)”)举行在[t,y] = ode45(@(t,y) 2*t， [0 3]， 0);情节(t y“o”)包含(“t”), ylabel (“y (t)”)标题(用ode45解y " =2t）

刚性微分方程

对于一些ODE问题，与积分区间相比，求解器所采取的步长被迫下降到一个不合理的小水平，即使在解曲线平滑的区域。这些步长可能非常小，以至于在短时间间隔内遍历可能需要进行数百万次计算。这可能会导致解算器的集成失败，但即使它成功了，也需要很长时间才能完成。

在ODE求解器中引起这种行为的方程被称为僵硬的．这是对一个事实的承认，即这些方程是顽固的，不容易用数值技术来计算。僵硬的ode带来的问题是显式求解器(例如数值)在达成解决方案方面的速度慢得令人难以忍受。这就是为什么数值被归类为一个非刚性求解器一起ode23而且ode113．这些求解器都很难对僵硬的方程进行积分。

方程刚度无法给出精确的定义，因为有几个因素会导致它。刚度由特定方程、所使用的ODE求解器、初始条件和求解器所使用的误差容限的组合产生。以下关于刚度的陈述，归因于Lambert[6]，在许多刚度ode的例子中得到了证明，但也存在反例，因此它们不是刚度的真实定义:

这些陈述的一个共同主题是，刚度可以由问题中某个地方的缩放差异引起。这种尺度上的差异(例如，如果雅可比矩阵$J = \partial f_n/\partial y_i$具有较大的负特征值比例)限制了求解器在执行积分时可以采取的步长。为了在解决方案中保持任何容错或稳定性的概念，微小的步长是必要的。

例如，描述化学反应的方程经常显示出刚度，因为溶液的组分通常在不同的时间尺度上变化很大(同时发生的反应既非常慢又非常快)。

但是，有一些求解器是专门设计来处理复杂的ode的。为刚性问题设计的求解器通常每一步要做更多的工作，结果是与非刚性求解器相比，它们能够采取更大的步骤，并享受更好的数值稳定性。刚性求解器是隐式的，因为计算$y$需要使用线性代数来求解线性方程组。雅可比矩阵用于估计ODE在积分过程中的局部行为，因此提供解析雅可比矩阵可以提高MATLAB的刚性ODE求解器的性能。

这只是对刚度的粗略处理，因为这是一个复杂的主题。看到常微分方程:刚度为了更深入地了解。

综上所述，MATLAB中的非刚性求解器为:

刚性求解器是(当数值是缓慢的):

应该注意的是，非刚性求解器是这样的工作在棘手的问题上，他们只是反应异常缓慢。类似地，为刚性问题设计的求解器也可以用于非刚性问题，但由于它们每一步要做更多的工作，因此在不必要的额外工作时，它们的效率比非刚性求解器要低。因此，方程刚度是求解器效率的问题，目标是在求解器在每一步中所做的工作和解的精度之间取得正确的平衡。

解决的建议

以下建议改编自MATLAB数学文档：

为该方法问题:

为僵硬的问题(数值是缓慢的):

下面的图表显示了基本的建议。在大多数情况下，在求解器中你需要做的唯一选择就是使用ode15s而不是数值．

例1:阻尼摆

阻尼摆的运动方程是，

$ $ \ ddot{\θ}= - \压裂{b} {m} \点{\θ}- \压裂{mg}{左(m-2b \右)L \} \罪\θ$ $

其中$g$是引力常数，$m$是弹子的质量，$L$是弦的长度，$b$是阻尼系数。我们的目标是求出摆偏离垂线的角度$\theta$，以及角度变化的速率$\theta '$。

一些自然的初始条件是$\theta_0 = \pi/4$和$\theta '_0 = 0$，这表明你在松开钟摆之前将钟摆举到45度角，并且它没有初始角速度。由于阻尼系数的存在，你可以预期钟摆会慢慢失去动量，然后回到静止状态。

该文件pendulumODE.m将问题重新表述为一阶ode的耦合系统:

数组$ $ \开始{}{cl} y_1’& = y_2 \ \ y_2 ' & = - \压裂{b} {m} y_2 - \压裂{mg} {L (m-2b)}罪(y_1) \{数组}$ $

然后用数值，ode15s，ode23,ode113．$y_1万博 尤文图斯 = \theta$的解被绘制出来，文件返回每个求解器的统计信息。在显示执行时间时，“显示的时间可能会有所不同”。

函数pendulumODE opts = odeset(“统计数据”，“上”）;Tspan = [0 25];Y0 = [pi/ 4,0];disp (＇ ＇), disp (ode45的统计数据:) tic， [t1,y1] = ode45(@pendode, tspan, y0, opts);toc disp (＇ ＇), disp (“ode15s的统计:”) tic， [t2,y2] = ode15s(@pendode, tspan, y0, opts);toc disp (＇ ＇), disp (ode23的统计数据:) tic， [t3,y3] = ode23(@pendode, tspan, y0, opts);toc disp (＇ ＇), disp (ode113的统计数据:) tic， [t4,y4] = ode113(@pendode, tspan, y0, opts);Toc图subplot(2,2,1)， plot(t1,y1(:，1)，“o”)， xlim([0 25])， title(“数值”) subplot(2,2,2)， plot(t2,y2(:，1)，“o”)， xlim([0 25])， title(“ode15s”) subplot(2,2,3)， plot(t3,y3(:，1)，“o”)， xlim([0 25])， title(“ode23”) subplot(2,2,4)， plot(t4,y4(:，1)，“o”)， xlim([0 25])， title(“ode113”）函数Dy2dt2 =摆管(t,y) g = 9.8;% m / s ^ 2M = 1;bob质量%L = 2;%钟摆长度，单位为米B = 0.2;%阻尼系数Dy2dt2 = [y(2);- b / m * y (2) - g / L *罪(y (1)));结束结束

pendulumODE

ode45的统计数据:75 successful steps 0 failed attempts 451 function evaluation运行时间为0.011970秒。ode15s的统计数据:183个成功步骤9个失败尝试315个函数计算1个偏导数31个LU分解311个线性系统的解耗时为0.081467秒。万博 尤文图斯ode23的统计数据:263个成功的步骤36个失败的尝试898个函数计算时间为0.015948秒。ode113的统计信息:159 successful steps 3 failed attempts 322 function evaluation运行时间为0.018056秒。

这些解都表现得很好，但阻尼摆是非刚性问题的一个很好的例子数值执行得很好。在这种情况下ode15s需要做额外的工作，以实现一个较差的解决方案。

例2:范德波尔振荡器

当非线性参数$\mu$较大时，范德堡尔振子方程在一定区间内变得刚性:

$$ ddot{y} - \mu \left(1-y^2\right)\dot{y}+y=0$$

这个方程的非线性完全包含在包含$\mu$的项中:注意，如果$\mu=0$，这个方程就退化为一个简单谐振子的方程，它具有规则的周期行为。

试图用数值遇到了严重的阻力，需要数百万次评估和30分钟以上的执行(我在35分钟后停止了执行)。由于问题很明显是刚性的，这个例子比较了刚性求解器。

该文件vanderpolODE.m找到$\mu=1000$的解使用ode15s，ode23s，ode23t,ode23tb．函数文件vdp1000.m并将该方程编码为一阶ode的耦合系统:

数组$ $ \开始{}{cl} y_1 ' & = y_2 \ \ y_2 ' & = \μ\离开(1-y_1 ^ 2 \右)y_2——y_1 \{数组}$ $

雅可比矩阵是用来辅助求解的，它的使用反映在偏导数计算的数量上。

函数vanderpolODE opts = odeset(“统计数据”，“上”，的雅可比矩阵, @J);Tspan = [0 3000];Y0 = [2 0];disp (＇ ＇), disp (“ode15s的统计:”) tic， [t1,y1] = ode15s(@vdp1000, tspan, y0, opts);toc disp (＇ ＇), disp (“ode23的统计数据:”) tic， [t2,y2] = ode23s(@vdp1000, tspan, y0, opts);toc disp (＇ ＇), disp (“ode23t的统计数据:”) tic， [t3,y3] = ode23t(@vdp1000, tspan, y0, opts);toc disp (＇ ＇), disp (ode23tb的统计数据:) tic， [t4,y4] = ode23tb(@vdp1000, tspan, y0, opts);Toc图subplot(2,2,1)， plot(t1,y1(:，1)，“o”)， ylim([-4 4])， title(“ode15s”) subplot(2,2,2)， plot(t2,y2(:，1)，“o”)、标题(“ode23s”) subplot(2,2,3)， plot(t3,y3(:，1)，“o”)、标题(“ode23t”) subplot(2,2,4)， plot(t4,y4(:，1)，“o”)、标题(“ode23tb”）函数dfdy = J(t,y) MU = 1000;非线性系数dfdy = [0 1 -2*MU*y(1)*y(2)-1 MU*(1-y(1)^2)];结束结束

vanderpolODE

ode15s的统计:591成功步骤225失败尝试1749函数计算45偏导数289 LU分解1747线性系统的解决方案耗时为0.192955秒。万博 尤文图斯ode23s的统计数据:741个成功步骤13个失败尝试2251个函数计算741个偏导数754个LU分解2262个线性系统的解耗时为0.092824秒。万博 尤文图斯ode23t的统计数据:776成功步骤94失败尝试2014函数评估36偏导数294 LU分解2012线性系统的解决方案运行时间为0.218996秒。万博 尤文图斯ode23tb的统计数据:573个成功步骤93个失败尝试2816个函数计算44个偏导数269个LU分解3415个线性系统的解耗时为0.202237秒。万博 尤文图斯

图是$y_1$的解。万博 尤文图斯对于这个问题，ode23s执行速度最快，失败的步骤最少。所提供的雅可比矩阵很有帮助ode23s求每一步的偏导数。ode23tb也用最少的步数解决了问题，表现优异ode15s．这个问题是一个很好的例子，一个生硬的问题与粗糙的公差ode23s而且ode23tb能表现得更好ode15s．但实际上，所有的硬求解器在这个问题上都表现得很好，并且相比之下节省了大量的时间数值．

总结

尽管所有的ODE求解器都能够处理相同的问题，但建议您从这里开始数值．然后，如果问题表现出僵硬的迹象，ode15s是不错的第二选择。然后，其他求解器根据特定问题的性质以及是否可以提供额外的信息(如雅可比矩阵)提供进一步的细化。

评论

你的工作涉及到MATLAB的ODE求解器的使用吗?如果有，分享你的经验在这里．

进一步的阅读

C.莫勒，常微分方程MATLAB数值计算电子版:The MathWorks, Inc.， Natick, MA, 2004

[2]香波，L. F.和M. W. Reichelt，MATLAB ODE套件SIAM科学计算杂志，Vol. 18, 1997。

C.莫勒，常微分方程:刚度克利夫角:克利夫·莫勒论数学与计算，2014

[4] l.f.香波，常微分方程的数值解查普曼和霍尔，1994年

香波，L. F.格拉德威尔，I.和S.汤普森，用MATLAB求解ode，剑桥大学出版社，2003年

j·d·兰伯特，常微分系统的数值方法，纽约:Wiley, 1992年

版权所有The MathWorks, Inc.

获取MATLAB代码

MATLAB®R2015a发布