解刚性常微分方程

打开脚本

这一页包含两个例子，解决刚性常微分方程使用ode15s．MATLAB®有四个为刚性ode设计的求解器。

ode15s
ode23s
ode23t
ode23tb

对于大多数棘手的问题，ode15s表现最好的。然而,ode23s，ode23t和ode23tb如果问题允许粗略的误差容限，则可以提高效率。

什么是僵硬ODE?

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

在ODE解算器中导致这种行为的方程称为僵硬的. 刚性ODE带来的问题是显式解算器（例如数值)在实现解决方案方面，他们的速度慢得无法忍受。这就是为什么数值被归类为非iff解算器随着奥德23，ode78，ode89和ode113．

为刚性常微分方程设计的解算器，称为刚性解算器，通常每一步做更多的工作。结果是，与非刚性求解器相比，它们能够采取更大的步骤，并提高了数值稳定性。

解算器选项

对于刚性问题，指定雅可比矩阵使用odeset特别重要。刚性解算器使用雅可比矩阵 $ p_partial f_i / p_partial y_j$ 在集成过程中估计ODE的局部行为，因此提供雅可比矩阵(或者，对于大型稀疏系统，提供其稀疏模式)对于效率和可靠性是至关重要的。使用雅可比矩阵，JPattern,或矢量化选择odeset来指定关于雅可比矩阵的信息。如果你不提供雅可比矩阵，那么求解器就会用有限差分数值估计它。

看到odeset有关其他解算器选项的完整列表。

例:Stiff van der Pol方程

范德波尔方程是一个二阶微分方程

$y''u 1-\mu\左（1-y''u 1^2\右）y''u 1+y''u 1=0$

在哪里 $美元\μ& # 62;0＄$ 是标量参数。当 $\μ= 1美元$ 由此产生的常微分方程组是非刚性的，并且易于使用数值. 然而，如果你增加 $\μ美元$ 到1000时，解发生了剧烈的变化，并在更长的时间尺度上表现出振荡。初值问题的近似解变得更加困难。因为这个特殊的问题是僵硬的，一个解决非僵硬问题的方法，例如数值它的效率太低，不实用。使用硬解算器，如ode15s而不是这个问题。

通过替换将范德波尔方程改写为一阶常微分方程组 $ y ' _1 = y_2美元 . 由此产生的一阶常微分方程系统为

$；\begin{array}{cl} ；y''u 1&；=y''u 2\\ ；y''u 2&；=\mu（1-y''u 1^2）y''u 2-y''u 1.\end{array}$

的vdp1000函数求范德波尔方程的值 $\μ= 1000美元$ ．

作用dydt = vdp1000 (t, y)计算mu = 1000时的范德堡尔ode。％%参见ODE15S, ODE23S, ODE23T, ODE23TB。亚采克·基尔赞卡和劳伦斯·f·沙宾版权所有1984-2014 The MathWorks, Inc.dydt = [y (2);1000 * (1 y (1) ^ 2) * y (2) - y (1)];

使用ode15s函数来解决初始条件为向量的问题[2; 0]，时间间隔为[0 3000]. 出于缩放原因，请仅绘制解决方案的第一个组件。

[t,y] = ode15s(@vdp1000，[0 3000]，[2;0]);情节(t y (: 1),“o”）;标题(van der Pol方程的解，\mu = 1000）;包含(“t”）;ylabel (“解决方案y_1”）;

的vdpode函数也解决了同样的问题，但它接受用户指定的值 $\μ美元$ ．方程变得越来越僵硬 $\μ美元$ 增加。

示例:稀疏布鲁塞尔系统

经典的Brusselator方程组可能是大的、僵硬的和稀疏的。Brusselator方程组模拟化学反应中的扩散，并由一个涉及美元$ ， $v$ ， $u'$ 和 $v'$ ．

$数组$ $ \开始{}{cl} u”_i & # 38; = 1 + u_i ^ 2 v_i-4u_i +α\ \离开(N + 1 \右)& # xA; ^ 2 \离开(u_张{}2 _i + u_ {i + 1} \) \ \ v ' _i & # 38; = 3 u_i-u_i ^ 2 v_i + \α# xA; \离开(N + 1 \右)^ 2 \离开(v_张{}- 2 v_i + v_ {i + 1} \) \{数组}$ $$

函数文件布鲁塞德在时间间隔上解这组方程[0,10]与 $\α= 1/50美元$ ．初始条件为

$\begin{array}{cl}u_j（0）和#38；=1+\sin（2\pix_j）\\v_j（0）和#38= ；3、\end{array}$

在哪里 x_j = i / N + 1美元为 i = 1美元,…,N ．因此,有 2 n美元而是雅可比矩阵 $ f / y$ 是否有一个宽度为5的带状矩阵，如果方程是有序的 u_1美元,v_1, u_2、v_2…美元像 N美元，问题变得越来越僵硬，雅可比矩阵变得越来越稀疏。

函数调用brussode (N),因为 $N\ge 2$ ，指定的值N在方程组中，对应于网格点的个数。默认情况下,布鲁塞德使用 $N=20$ ．

布鲁塞德包含几个子功能：

嵌套函数f（t，y）对布鲁塞尔振子问题的方程组进行编码，返回一个向量。
当地的函数jpattern (N)返回一个1和0的稀疏矩阵，显示非零在雅可比矩阵中的位置。此矩阵被指定给JPattern选项结构的字段。ODE解算器使用此稀疏模式以稀疏矩阵的形式生成雅可比矩阵。在问题中提供这种稀疏模式可以显著减少生成2N×2N雅可比矩阵所需的函数求值数量，从2N个求值减少到4个。

作用brussode (N)模拟化学反应的困难问题。%参数N >= 2用于指定网格点的个数;的生成的系统由2N个方程组成。缺省情况下，N为20。的%随着N的增加，问题变得越来越僵硬，越来越稀疏%增加。这个问题的雅可比矩阵是一个稀疏常数矩阵%(带状带宽为5)。％%属性“JPattern”用于为解算器提供稀疏由1和0组成的%矩阵，表示雅可比矩阵中非零的位置%df/dy。默认情况下，ODE套件的刚性解算器生成雅可比矩阵%以全矩阵的形式进行数值计算。但是，当稀疏模式%，求解器使用它来生成雅可比矩阵的数值形式为a%稀疏矩阵。提供稀疏模式可以显著减少生成雅可比矩阵所需的函数值的百分比%加速整合。对于BRUSSODE问题，只有4个评价%函数需要计算2N x 2N雅可比矩阵。％%设置“矢量化”属性表示函数f%矢量化。％% E.海勒，G.万纳，求解常微分方程II，僵硬和微分代数问题，施普林格-弗拉格，柏林，%1991年，第5-8页。％%参见ODE15S, ODE23S, ODE23T, ODE23TB, ODESET, FUNCTION_HANDLE。%马克·W·雷切尔和劳伦斯·F·萨姆平，1994年8月30日版权所有1984-2014 The MathWorks, Inc.问题参数，与嵌套函数共享。如果nargin<1 N = 20;结束tspan = [0;10);罪y0 =(1 +(2 *π/ (N + 1) * (1: N));repmat (1, N)];选择= odeset (矢量化的，“上”，“JPattern”，jpattern（N））；[t，y]=ode15s（@f，tspan，y0，options）；u=y（：，1:2:end）；x=（1:N）/（N+1）；图；冲浪（x，t，u）；视图（-40,30）；xlabel(“空间”）;ylabel (“时间”); 兹拉贝尔(解决你的); 头衔([“为N的布鲁塞尔人”num2str（N）]；% -------------------------------------------------------------------------%嵌套函数——N由外部函数提供。％作用dydt = f (t, y)%导数函数c=0.02*（N+1）^2；dydt=0（2*N，大小（y，2））；%预分配dy/dt%计算网格一侧的函数的两个分量%(带边缘条件)。i=1；dydt（i，：）=1+y（i+1，：）.*y（i，：）。^2-4*y（i，：）+c*（1-2*y（i，：）+y（i+2，：）；dydt（i+1，：）=3*y（i，：）-y（i+1，：）.*y（i，：）。^2+c*（3-2*y（i+1，：）+y，：）；在所有内部网格点上计算函数的2个分量。我= 3:2:2 * n - 3;= 1 + y(i+1，:).*y(i，:)。^2 - 4*y(i，:) +．..c*（y（i-2，：）-2*y（i，：）+y（i+2，：）；dydt（i+1，：）=3*y（i，：）-y（i+1，：）.*y（i，：）。^2+．..: c * (y(张)2 *(我+ 1,)+ y (i + 3,:));%在网格的另一边评估函数的两个组件%(带边缘条件)。i=2*N-1；dydt（i，：）=1+y（i+1，：）.*y（i，：）。^2-4*y（i，：）+c*（y（i-2，：）-2*y（i，：）+1）；dydt（i+1，：）=3*y（i，：）-y（i+1，：）.*y（i，：）.^2+c*（y（i-1，：）-2*y（i+1，：）+3）；结束% -------------------------------------------------------------------------结束%布鲁塞德% ---------------------------------------------------------------------------%子函数——稀疏模式％作用S = jpattern (N)%雅可比稀疏模式B=一（2*N，5）；B（2:2:2*N，2）=零（N，1）；B（1:2:2*N-1,4）=零（N，1）；S=spdiags（B，-2:2,2*N，2*N）；结束% ---------------------------------------------------------------------------