解算刚性 ODE
本页包含两个使用 ode15s
解算刚性常微分方程的示例。MATLAB® 拥有四个专用于刚性 ODE 的求解器。
ode15s
ode23s
ode23t
ode23tb
对于大多数刚性问题,ode15s
的性能最佳。但如果问题允许较宽松的误差容限,则 ode23s
、ode23t
和 ode23tb
可能更加高效。
什么是刚性 ODE?
对于一些 ODE 问题,求解器采用的步长被强制缩小为与积分区间相比过小的级别,甚至在解曲线平滑的区域亦如此。这些步长可能过小,以至于遍历很短的时间区间都可能需要数百万次计算。这可能导致求解器积分失败,即使积分成功也需要花费很长时间。
导致 ODE 求解器出现此行为的方程称为刚性方程。刚性 ODE 造成的问题是,显式求解器(例如 ode45
)获取解的速度慢得令人无法忍受。这是将 ode45
与 ode23
、ode78
、ode89
和 ode113
一同归类为非刚性求解器的原因所在。
专用于刚性 ODE 的求解器称为刚性求解器,它们通常在每一步中完成更多的计算工作。这样做的好处是,它们能够采用大得多的步长,并且与非刚性求解器相比提高了数值稳定性。
求解器选项
对于刚性问题,使用 odeset
指定雅可比矩阵尤为重要。刚性求解器使用雅可比矩阵 来预测 ODE 在积分过程中的局部行为,因此提供雅可比矩阵(或者对于大型稀疏方程组提供其稀疏模式)对于提高效率和可靠性而言至关重要。使用
odeset
的 Jacobian
、JPattern
或 Vectorized
选项来指定雅可比的相关信息。如果没有提供雅可比,则求解器将使用有限差分对其进行数值预测。
有关其他求解器选项的完整列表,请参阅 odeset
。
示例:刚性 van der Pol 方程
van der Pol 方程为二阶 ODE
其中 为标量参数。当
时,生成的 ODE 方程组为非刚性方程组,可以使用
ode45
轻松求解。但如果将 增大至 1000,则解会发生显著变化,并会在明显更长的时间段中显示振荡。求初始值问题的近似解变得更加复杂。由于此特定问题是刚性问题,因此专用于非刚性问题的求解器(如
ode45
)的效率非常低下且不切实际。针对此问题应改用 ode15s
等刚性求解器。
通过执行代换 ,将该 van der Pol 方程重写为一阶 ODE 方程组。生成的一阶 ODE 方程组为
vdp1000
函数使用 计算 van der Pol 方程。
function dydt = vdp1000(t,y) %VDP1000 Evaluate the van der Pol ODEs for mu = 1000. % % See also ODE15S, ODE23S, ODE23T, ODE23TB. % Jacek Kierzenka and Lawrence F. Shampine % Copyright 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]); plot(t,y(:,1),'-o'); title('Solution of van der Pol Equation, \mu = 1000'); xlabel('Time t'); ylabel('Solution y_1');
vdpode
函数也可以求解同一问题,但它接受的是用户指定的 值。随着
的增大,该方程组的刚性逐渐增强。
示例:稀疏 Brusselator 方程组
经典 Brusselator 方程组可能为大型刚性稀疏矩阵。Brusselator 方程组可模拟化学反应中的扩算,并表示为涉及 、
、
和
的方程组。
函数文件 brussode
使用 在时间区间
[0,10]
上对这组方程进行求解。初始条件为
其中,当 时,
。因此,该方程组中有
个方程,但如果这些方程按
的形式排序,则雅可比
为具有常量宽度 5 的带状矩阵。随着
的增大,此问题的刚性逐渐增强,并且雅可比矩阵的稀疏性也逐渐增大。
函数调用 brussode(N)
(其中 )为方程组中的
N
(对应于网格点数量)指定值。默认情况下,brussode
使用 。
brussode
包含一些子函数:
嵌套函数
f(t,y)
用于编写 Brusselator 问题的方程组代码,并返回一个向量。局部函数
jpattern(N)
返回由 1 和 0 组成的稀疏矩阵,从而显示雅可比矩阵中非零值的位置。此矩阵将赋给 options 结构体的JPattern
字段。ODE 求解器使用此稀疏模式,生成稀疏矩阵形式的雅可比数值矩阵。在问题中提供此稀疏模式可将生成 2N×2N 雅可比矩阵所需的函数计算量从 2N 次大幅减少至仅仅 4 次。
function brussode(N) %BRUSSODE Stiff problem modelling a chemical reaction (the Brusselator). % The parameter N >= 2 is used to specify the number of grid points; the % resulting system consists of 2N equations. By default, N is 20. The % problem becomes increasingly stiff and increasingly sparse as N is % increased. The Jacobian for this problem is a sparse constant matrix % (banded with bandwidth 5). % % The property 'JPattern' is used to provide the solver with a sparse % matrix of 1's and 0's showing the locations of nonzeros in the Jacobian % df/dy. By default, the stiff solvers of the ODE Suite generate Jacobians % numerically as full matrices. However, when a sparsity pattern is % provided, the solver uses it to generate the Jacobian numerically as a % sparse matrix. Providing a sparsity pattern can significantly reduce the % number of function evaluations required to generate the Jacobian and can % accelerate integration. For the BRUSSODE problem, only 4 evaluations of % the function are needed to compute the 2N x 2N Jacobian matrix. % % Setting the 'Vectorized' property indicates the function f is % vectorized. % % E. Hairer and G. Wanner, Solving Ordinary Differential Equations II, % Stiff and Differential-Algebraic Problems, Springer-Verlag, Berlin, % 1991, pp. 5-8. % % See also ODE15S, ODE23S, ODE23T, ODE23TB, ODESET, FUNCTION_HANDLE. % Mark W. Reichelt and Lawrence F. Shampine, 8-30-94 % Copyright 1984-2014 The MathWorks, Inc. % Problem parameter, shared with the nested function. if nargin<1 N = 20; end tspan = [0; 10]; y0 = [1+sin((2*pi/(N+1))*(1:N)); repmat(3,1,N)]; options = odeset('Vectorized','on','JPattern',jpattern(N)); [t,y] = ode15s(@f,tspan,y0,options); u = y(:,1:2:end); x = (1:N)/(N+1); figure; surf(x,t,u); view(-40,30); xlabel('space'); ylabel('time'); zlabel('solution u'); title(['The Brusselator for N = ' num2str(N)]); % ------------------------------------------------------------------------- % Nested function -- N is provided by the outer function. % function dydt = f(t,y) % Derivative function c = 0.02 * (N+1)^2; dydt = zeros(2*N,size(y,2)); % preallocate dy/dt % Evaluate the 2 components of the function at one edge of the grid % (with edge conditions). 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(i+3,:)); % Evaluate the 2 components of the function at all interior grid points. i = 3:2:2*N-3; dydt(i,:) = 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(i-1,:)-2*y(i+1,:)+y(i+3,:)); % Evaluate the 2 components of the function at the other edge of the grid % (with edge conditions). 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); end % ------------------------------------------------------------------------- end % brussode % --------------------------------------------------------------------------- % Subfunction -- the sparsity pattern % function S = jpattern(N) % Jacobian sparsity pattern B = ones(2*N,5); B(2:2:2*N,2) = zeros(N,1); B(1:2:2*N-1,4) = zeros(N,1); S = spdiags(B,-2:2,2*N,2*N); end % ---------------------------------------------------------------------------
通过运行函数 brussode
,对 时的 Brusselator 方程组求解。
brussode
通过为 brussode
指定输入,对 时的方程组求解。
brussode(50)
另请参阅
ode15s
| ode23s
| ode23t
| ode23tb