Main Content

本页的翻译已过时。点击此处可查看最新英文版本。

求解具有多个初始条件的 ODE 方程组

此示例比较求解具有多组初始条件的常微分方程组的两种方法。这些方法是:

  • 使用 for 循环执行多次仿真,每组初始条件对应一次仿真。此方法使用起来很简单,但对于大型方程组不能实现最优性能。

  • 向量化 ODE 函数,以同时针对各组初始条件求解方程组。对于大型方程组来说,这是更快的方法,但需要重写 ODE 函数以便正确地重构输入。

用于说明这些方法的方程是众所周知的 Lotka-Volterra 方程,它们是描述捕食者和猎物种群的一阶非线性微分方程。

问题描述

Lotka-Volterra 方程是由两个一阶非线性 ODE 组成的方程组,用于描述生物系统中捕食者和猎物的种群。随着时间的推移,捕食者和猎物的数量会根据方程发生变化

dxdt=αx-βxy,dydt=δxy-γy.

这些方程中的变量包括

  • x 是猎物的种群大小

  • y 是捕食者的种群大小

  • t 是时间

  • αβδγ 是描述两个物种之间交互的常量参数。此示例使用参数值 α=γ=1β=0.01δ=0.02

对于此问题,xy 的初始值是初始种群大小。求解这些方程,就能提供关于物种交互时种群如何随时间变化的信息。

求解具有一个初始条件的方程

要在 MATLAB 中求解 Lotka-Volterra 方程,请编写一个函数来对方程进行编码,指定积分的时间区间,并指定初始条件。然后,您可以使用一个 ODE 求解器(如 ode45)对方程组随时间的变化进行仿真。

对方程进行编码的函数是

function dpdt = lotkaODE(t,p)
% LOTKA Lotka-Volterra predator-prey model
delta = 0.02;
beta = 0.01;

dpdt = [p(1) .* (1 - beta*p(2));
        p(2) .* (-1 + delta*p(1))];
end

(该函数作为局部函数包含在示例的末尾。)

由于方程组中有两个方程,dpdt 是向量,其中每个元素对应一个方程。此外,对于每个解分量,解向量 p 中都有一个对应的元素:p(1) 表示原始方程中的 xp(2) 表示原始方程中的 y

接下来,将积分的时间区间指定为 [0,15],并将 xy 的初始种群大小设置为 50。

t0 = 0;
tfinal = 15;
p0 = [50; 50];

通过指定 ODE 函数、时间跨度和初始条件,用 ode45 求解方程组。绘制结果种群对时间的图。

[t,p] = ode45(@lotkaODE,[t0 tfinal],p0);
plot(t,p)
title('Predator/Prey Populations Over Time')
xlabel('t')
ylabel('Population')
legend('Prey','Predators')

Figure contains an axes. The axes with title Predator/Prey Populations Over Time contains 2 objects of type line. These objects represent Prey, Predators.

由于解呈现周期性,在相位图中绘制各解之间相对关系的图。

plot(p(:,1),p(:,2))
title('Phase Plot of Predator/Prey Populations')
xlabel('Prey')
ylabel('Predators')

Figure contains an axes. The axes with title Phase Plot of Predator/Prey Populations contains an object of type line.

所得的图显示针对给定初始种群大小的解。要求解针对不同初始种群大小的方程,请更改 p0 中的值,然后重新运行仿真。然而,此方法一次只能针对一个初始条件求解方程。接下来的两节说明针对许多不同初始条件求解的方法。

方法 1:用 for- 循环计算多个初始条件

求解具有多个初始条件的 ODE 方程组的最简单方法是使用 for 循环。此方法使用与单一初始条件方法相同的 ODE 函数,但 for 循环会自动执行求解过程。

例如,您可以将 x 的初始种群大小保持为 50,并使用 for 循环使 y 的初始种群大小在 10 到 400 之间变化。为 y0 创建一个种群大小向量,然后使用这些值执行循环以针对每组初始条件求解方程。用所有迭代的结果绘制一个相位图。

y0 = 10:10:400;
for k = 1:length(y0)
    [t,p] = ode45(@lotkaODE,[t0 tfinal],[50 y0(k)]);
    plot(p(:,1),p(:,2))
    hold on
end
title('Phase Plot of Predator/Prey Populations')
xlabel('Prey')
ylabel('Predators')
hold off

Figure contains an axes. The axes with title Phase Plot of Predator/Prey Populations contains 40 objects of type line.

该相位图显示不同初始条件组的所有计算结果。

方法 2:用向量化 ODE 函数计算多个初始条件

求解多初始条件的 ODE 方程组的另一种方法是重写 ODE 函数,以便同时求解所有方程。其实现步骤如下:

  • 将所有初始条件作为一个矩阵提供给 ode45。该矩阵的大小为 s×n,其中 s 是解分量数目,n 是要求解的方程的初始条件的数目。矩阵中的每列表示方程组的一组完整初始条件。

  • ODE 函数必须接受额外输入参数 n,即初始条件的数目。

  • 在 ODE 函数中,求解器将解分量 p 作为列向量传递。ODE 函数必须将向量重构为 s×n 大小的矩阵。这样,矩阵的每行都包含每个变量的所有初始条件。

  • ODE 函数必须求解向量化格式的方程,以便表达式接受向量作为解分量。换句话说,f(t,[y1 y2 y3 ...]) 必须返回 [f(t,y1) f(t,y2) f(t,y3) ...]

  • 最后,ODE 函数必须将其输出重构回向量,以便 ODE 求解器从每次函数调用接收一个向量。

如果您按照这些步骤操作,则 ODE 求解器可以使用解分量的向量来求解方程组,而 ODE 函数将向量重构为矩阵,并针对所有初始条件求解每个解分量。结果是,您可以在一次仿真中针对所有初始条件求解方程组。

要对 Lotka-Volterra 方程组实现此方法,首先要找到初始条件的数量 n,然后构成一个初始条件矩阵。

n = length(y0);
p0_all = [50*ones(n,1) y0(:)]';

接下来,重写 ODE 函数,使其接受 n 作为输入。使用 n 将解向量重构为矩阵,然后求解向量化方程组,并将输出重构回向量。经修改后执行这些任务的 ODE 函数是

function dpdt = lotkasystem(t,p,n)
%LOTKA  Lotka-Volterra predator-prey model for system of inputs p.
delta = 0.02;
beta = 0.01;

% Change the size of p to be: Number of equations-by-number of initial
% conditions.
p = reshape(p,[],n);

% Write equations in vectorized form.
dpdt = [p(1,:) .* (1 - beta*p(2,:)); 
        p(2,:) .* (-1 + delta*p(1,:))];

% Linearize output.
dpdt = dpdt(:);
end

使用 ode45 针对所有初始条件求解方程组。由于 ode45 要求 ODE 函数接受两个输入,请使用匿名函数将 n 的值从工作区传入 lotkasystem

[t,p] = ode45(@(t,p) lotkasystem(t,p,n),[t0 tfinal],p0_all);

将输出向量重构为 (numTimeSteps*s)×n 大小的矩阵。输出 p(:,k) 的每列都包含针对一组初始条件的解。绘制解分量的相位图。

p = reshape(p,[],n);
nt = length(t);
for k = 1:n
    plot(p(1:nt,k),p(nt+1:end,k))
    hold on
end
title('Predator/Prey Populations Over Time')
xlabel('t')
ylabel('Population')
hold off

Figure contains an axes. The axes with title Predator/Prey Populations Over Time contains 40 objects of type line.

结果与通过 for 循环方法获得的值类似。但是,您应牢记向量化解方法的一些特性:

  • 计算的解可能与从单个初始输入计算的解略有不同。出现差异是因为 ODE 求解器对整个方程组应用范数检查来计算时间步的大小,因此解的时间步进行为略有不同。时间步的变化通常不会影响解的精确度,而是影响计算解的时间点。

  • 对于自动计算方程组的数值 Jacobian 矩阵的刚性 ODE 求解器(ode15sode23sode23tode23tb),使用 odesetJPattern 选项指定 Jacobian 矩阵的分块对角稀疏模式可以提高计算效率。分块对角 Jacobian 矩阵是重写的 ODE 函数对输入重构后产生的。

比较计时结果

使用 timeit 对上述每种方法进行计时。求解含一组初始条件的方程的计时结果将作为基线值来衡量各方法。

% Time one IC
baseline = timeit(@() ode45(@lotkaODE,[t0 tfinal],p0),2);

% Time for-loop
for k = 1:length(y0)
    loop_timing(k) = timeit(@() ode45(@lotkaODE,[t0 tfinal],[50 y0(k)]),2);
end
loop_timing = sum(loop_timing);

% Time vectorized fcn
vectorized_timing = timeit(@() ode45(@(t,p) lotkasystem(t,p,n),[t0 tfinal],p0_all),2);

创建一个计时结果表。将所有结果乘以 1e3,以毫秒为单位表示时间。其中一列为每个解的平均用时,该值通过将计时时间除以初始条件数目求得。

TimingTable = table(1e3.*[baseline; loop_timing; vectorized_timing], 1e3.*[baseline; loop_timing/n; vectorized_timing/n],...
    'VariableNames',{'TotalTime (ms)','TimePerSolution (ms)'},'RowNames',{'One IC','Multi ICs: For-loop', 'Mult ICs: Vectorized'})
TimingTable=3×2 table
                            TotalTime (ms)    TimePerSolution (ms)
                            ______________    ____________________

    One IC                      1.4748               1.4748       
    Multi ICs: For-loop         34.158              0.85395       
    Mult ICs: Vectorized        6.8257              0.17064       

TimePerSolution 列显示向量化方法是三种方法中最快的。

局部函数

此处列出了 ode45 为计算解而调用的局部函数。

function dpdt = lotkaODE(t,p)
% LOTKA Lotka-Volterra predator-prey model
delta = 0.02;
beta = 0.01;

dpdt = [p(1) .* (1 - beta*p(2));
        p(2) .* (-1 + delta*p(1))];
end
%------------------------------------------------------------------
function dpdt = lotkasystem(t,p,n)
%LOTKA  Lotka-Volterra predator-prey model for system of inputs p.
delta = 0.02;
beta = 0.01;

% Change the size of p to be: Number of equations-by-number of initial
% conditions.
p = reshape(p,[],n);

% Write equations in vectorized form.
dpdt = [p(1,:) .* (1 - beta*p(2,:)); 
        p(2,:) .* (-1 + delta*p(1,:))];

% Linearize output.
dpdt = dpdt(:);
end

另请参阅

|

相关主题