求解具有多初始条件的ode系统

打开生活的脚本

这个例子比较了两种方法来求解具有多组初始条件的常微分方程组。的技术是:

使用一个为-loop执行几个模拟，每个初始条件集一个。这种技术使用简单，但不能为大型系统提供最佳性能。
对ODE函数进行矢量化，同时求解所有初始条件集的方程组。这种技术对于大型系统来说是更快的方法，但是需要重写ODE函数，以便它能够正确地重新定义输入。

用来证明这些技术的方程是著名的Lotka-Volterra方程，这是一阶非线性微分方程，描述了捕食者和猎物的数量。

问题描述

的生态方程是一个由两个一阶非线性ode组成的系统，用于描述生物系统中捕食者和猎物的数量。随着时间的推移，捕食者和猎物的数量会根据这个方程发生变化

$\begin{array}{l} \frac{dx}{dt} ＝ α x - β xy ， \\ \frac{dy}{dt} ＝ δ xy - γ y ． \end{array}$

这些方程中的变量是

$x$ 猎物的数量是多少
$y$ 捕食者的数量是多少
$t$ 是时候
$α$ ， $β$ ， $δ$ , $γ$ 是描述两种物种间相互作用的常数参数。本示例使用参数值 $α ＝ γ ＝ 1$ ， $β ＝ 0 ． 01$ , $δ ＝ 0 ． 02$ ．

对于这个问题，初始值 $x$ 和 $y$ 为初始种群大小。解出这个方程，就可以了解种群如何随着时间的推移而随着物种的相互作用而变化。

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

要在MATLAB中求解Lotka-Volterra方程，需要编写一个函数对方程进行编码，指定积分的时间间隔，并指定初始条件。然后您可以使用一个ODE求解器，例如数值，以模拟系统随时间的变化。

编码方程的函数是

函数dpdt = lotkaODE (t, p)% LOTKA LOTKA - volterra捕食-猎物模型δ= 0.02;β= 0.01;DPDT = [p(1) .* (1 - beta*p(2))];P (2) .* (-1 + delta* P (1))];结束

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

因为系统中有两个方程，dpdt是每个方程都有一个元素的向量。还有解向量p每个解决方案组件都有一个元素:(1页)代表 $x$ 在原始方程中，和(2页)代表 $y$ 在原始方程中。

接下来，指定积分的时间间隔为 $［ 0 ， 15 ］$ 并设置初始种群大小 $x$ 和 $y$ 到50。

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

用下列方法解方程组数值通过指定ODE函数、时间跨度和初始条件。绘制结果种群与时间的关系图。

[t,p] = ode45(@lotkaODE，[t0 tfinal]，p0);情节(t, p)标题(“随着时间的推移，捕食者/猎物的数量”)包含(“t”) ylabel (“人口”)传说(“猎物”，“捕食者”）

图中包含一个轴对象。标题为Predator/Prey population Over Time的轴对象包含2个类型为line的对象。这些物体代表猎物，捕食者。

由于解具有周期性，所金宝搏官方网站以在相图中绘制解与解的相对关系。

情节(p(: 1),(2):,)标题(“捕食者/猎物种群的阶段图”)包含(“猎物”) ylabel (“捕食者”）

图中包含一个轴对象。标题为“捕食者/猎物种群相图”的轴对象包含一个类型为线的对象。

结果图显示了给定初始种群大小的解。为了求解不同初始种群大小的方程组，改变中的值p0然后重新运行模拟。然而，这种方法一次只解一个初始条件下的方程组。接下来的两部分将描述解决许多不同初始条件的技术。

方法一:用`对- - - - - -`循环

求解具有多个初始条件的ode方程组最简单的方法是用为循环。这种技术使用与单一初始条件技术相同的ODE函数，但是为-loop使解决过程自动化。

例如，可以保存初始种群大小 $x$ 常数为50，并使用为-loop来改变初始种群大小 $y$ 在10到400之间。为创建一个人口大小矢量y0，然后循环遍历这些值来求解每一组初始条件的方程。用所有迭代的结果绘制一个阶段图。

y0 = 10:10:400;为k = 1:长度(y0) (t, p) =数值(@lotkaODE [t0 tfinal], [50 y0 (k)]);情节(p (: 1), (2):,)在结束标题(“捕食者/猎物种群的阶段图”)包含(“猎物”) ylabel (“捕食者”)举行从

图中包含一个轴对象。标题为“捕食者/猎物种群相图”的轴对象包含40个直线类型的对象。

相图显示了不同初始条件集的所有计算解。金宝搏官方网站

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

求解具有多个初始条件的ODE方程组的另一种方法是重写ODE函数，使所有方程组同时求解。做到这一点的步骤是:

提供所有的初始条件数值作为一个矩阵。矩阵的大小是年代——- - - - - -n,在那里年代解决方案组件和的数量是多少n是初始条件的个数。矩阵中的每一列代表系统的一组完整的初始条件。
ODE函数必须接受额外的输入参数n，初始条件的个数。
在ODE函数内部，求解器将解组件p作为列向量传递。ODE函数必须将向量重塑为具有大小的矩阵年代×n。矩阵的每一行包含了每个变量的所有初始条件。
ODE函数必须以矢量格式求解方程，以便表达式接受向量作为解的组成部分。换句话说,F (t，[y1 y2 y3…])必须返回[f(t,y1) f(t,y2) f(t,y3)…]．
最后，ODE函数必须将其输出重新设置为向量，以便ODE求解器从每次函数调用中接收回一个向量。

如果您遵循这些步骤，那么ODE解算器可以使用一个向量作为解的分量来解方程组，而ODE函数将向量重塑为一个矩阵，并解决所有初始条件下的每个解的分量。结果是你可以在一个模拟中解决系统的所有初始条件。

要实现Lotka-Volterra系统的这种方法，首先要找到初始条件的数量n，然后形成一个初始条件矩阵。

n =长度(y0);P0_all = [50*ones(n,1) y0(:)]';

接下来，重写ODE函数，使其接受n作为输入。使用n将解向量重塑为矩阵，然后求解矢量化系统并将输出重塑为向量。执行这些任务的修改过的ODE函数是

函数dpdt = lotkasystem (t, p, n)输入系统的LOTKA - volterra捕食-食饵模型。δ= 0.02;β= 0.01;%改变p的大小为%的条件。[] p =重塑(p, n);把方程写成向量化的形式。DPDT = [p(1，:) .* (1 - beta*p(2，:));P (2，:) .* (-1 + * P (1，:))];%线性化输出。dpdt = dpdt (:);结束

解方程组的所有初始条件使用数值．自数值要求ODE函数接受两个输入，使用匿名函数传递值n从工作区到lotkasystem．

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

将输出向量重塑为具有大小的矩阵(numTimeSteps *年代)——- - - - - -n．输出的每一列p (:, k)包含一组初始条件的解。金宝搏官方网站绘制溶液组分的相图。

[] p =重塑(p, n);元=长度(t);为K = 1:n plot(p(1:nt, K)，p(nt+1:end, K)) hold在结束标题(“随着时间的推移，捕食者/猎物的数量”)包含(“t”) ylabel (“人口”)举行从

图中包含一个轴对象。标题为Predator/Prey Populations Over Time的轴对象包含40个line类型的对象。

所得结果与该方法所得结果具有可比性为循环技术。然而，你应该记住矢量化解法的一些特性:

计算出的解可能与从单个初始输入金宝搏官方网站计算出的解略有不同。出现这种差异是因为ODE求解器对整个系统应用规范检查来计算时间步长，所以解的时间步长行为略有不同。时间步长的变化通常不会影响解的精度，而是影响解的计算时间。
对于刚性ODE求解器(ode15s，ode23s，ode23t，ode23tb)，自动计算系统的数值雅可比矩阵，指定雅可比矩阵的块对角稀疏模式JPattern选择odeset可以提高计算效率。雅可比矩阵的块对角形式来自于在重写的ODE函数中执行的输入整形。

比较时间结果

计算前面每个方法的使用时间时间．用一组初始条件解方程组的时间作为基线数，以查看方法的规模。

% Time one IC基线= timeit(@() ode45(@ lotkaode，[t0 tfinal]，p0)，2);%的时间循环为@() ode45(@ lotkaode，[t0 tfinal]，[50 y0(k)])，2);结束loop_timing =总和(loop_timing);%时间矢量化fcnVectorized_timing = timeit(@() ode45(@(t,p) lotkasystem(t,p,n)，[t0 tfinal]，p0_all)，2);

创建一个包含计时结果的表。将所有结果乘以1e3以毫秒表示时间。包括一个列，列中包含每个解所需的时间，这个时间除以要解的初始条件的数量。

TimingTable =表(1 e3。*(基线;loop_timing;vectorized_timing], 1 e3。*[基线;loop_timing / n;vectorized_timing / n),．..“VariableNames”, {“TotalTime (ms)”，“TimePerSolution (ms)”}，“RowNames”, {一个集成电路的，“多ICs:循环”，“乘ICs:矢量化”}）

TimingTable =3×2表TotalTime (ms) TimePerSolution (ms) ______________ ____________________ One IC 1.6558 1.6558 Multi ICs: For-loop 34.333 0.5832 Multi ICs: Vectorized 5.1357 0.12839 . png (13.688 kb，下载次数:0

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

本地函数

这里列出的是本地函数数值调用来计算解决方案。金宝搏官方网站

函数dpdt = lotkaODE (t, p)% LOTKA LOTKA - volterra捕食-猎物模型δ= 0.02;β= 0.01;DPDT = [p(1) .* (1 - beta*p(2))];P (2) .* (-1 + delta* P (1))];结束%------------------------------------------------------------------函数dpdt = lotkasystem (t, p, n)输入系统的LOTKA - volterra捕食-食饵模型。δ= 0.02;β= 0.01;%改变p的大小为%的条件。[] p =重塑(p, n);把方程写成向量化的形式。DPDT = [p(1，:) .* (1 - beta*p(2，:));P (2，:) .* (-1 + * P (1，:))];%线性化输出。dpdt = dpdt (:);结束