基于常微分方程初值问题的ModelCoder固定步长求解器分析

ModelCoder解算器库中主要有两种类型的解算器——固定步长解算器和可变步长解算器。在这里主要讨论固定步长解算器。显而易见的，所谓固定步长解算器，就是步长大小在仿真过程中保持不变。

固定步长离散解算器通过将当前时间加上固定步长大小来计算下一个仿真步的时间，执行这种计算时，它们依赖于模型中的每个模块来更新其各个离散状态。它们不计算连续状态。仿真的准确性和时间长短取决于仿真执行的步长大小：步长越小，结果越准确，但仿真时间越长。默认情况下，ModelCoder会选择步长大小，但用户也可以自己选择步长大小。如果选择默认设置 “auto”，而且模型具有离散采样时间，则 ModelCoder会将步长大小设置为等于模型的基础采样时间。否则，如果不存在离散速率，ModelCoder会将步长大小设置为等于仿真开始时间与结束时间之差再除以 50。与固定步长离散解算器一样，固定步长连续解算器也通过将当前时间加上固定大小的时间步来计算下一个仿真时间。对于每一步，连续解算器使用数值积分方法来计算模型的连续状态值。计算这些值需要使用上一个时间步的连续状态，以及当前时间步与上一个时间步之间的中间点（子时间步）的状态导数。因此，固定步长连续解算器可以处理同时包含连续状态和离散状态的模型。本文将对固定步长连续解算器展开分析。

固定步长连续解算器使用数值积分方法通常涉及到常微分方程数值求解，一般来说为求解常微分方程初值问题，即寻找函数 y(t) 使之满足如下常微分方程：

所谓数值解法，就是求y(x)在离散结点xn处的函数近似值yn的方法，yn≈y(xn)。这些近似值称为常微分方程初值问题的数值解。相邻两个结点之间的距离称为步长。微分方程数值解的基本思想是：求解区间和方程离散化。

只要f(x,y)在[a,b]*R^d上连续，且关于y满足Lipschuitz条件，即存在与x,y无关的常数L使 |f(x,y1)-f(x,y2)|≤L|y1-y2| 对任意定义在[a,b]上的y1(x)和y2(x)都成立，则上述IVP(Initial-Value Problem)存在唯一解。

要计算出解函数y(x)在一系列节点 a=x0<x1<...<xn=b 处的近似值yn≈y(xn),n=1,2,...,N，其中y0=y(a)已知，节点间距位步长，通常采取等距网格，可以得到：

接着需要将微分方程离散化，在这里通过欧拉方法来举例介绍，通过向前差商近似导数得到：

移项变化可得：

将y0+hf(x0,y0)记为y1，可以得到一般的递推公式：

下面在ModelCoder中通过两个实例来解释说明，首先是将常数作为输入信号连接一个积分器，如图1所示，ModelCoder将 Integrator 模块作为具有一种状态的动态系统进行处理。模块动态由以下方程指定：

其中，u是模块输入，y是模块输出，x是模块状态，x0

是x的初始条件。

▲图1 积分器模型

虽然这些方程定义了连续时间下的准确关系，但ModelCoder使用数值逼近方法以有限精度来进行计算。ModelCoder可以使用若干不同的数值积分方法来计算模块的输出，每种方法都在特定的应用中各具优势。本例使用欧拉方法来计算模块的输出，具体地，在设置对话框的求解器窗格中可以选择欧拉方法所选求解器会使用当前输入值和前一个时间步的状态值计算 Integrator模块在当前时间步的输出。为支持此计算模型，Integrator模块会保存在当前时间步的输出，以供求解器计算其在下一个时间步的输出。

在设置中将固定步长设置“auto”，而且本例模型不存在离散速率，ModelCoder会将步长大小设置为等于仿真开始时间与结束时间之差再除以50，即10/50。

ModelCoder固定步长连续显式求解器通过状态和状态导数的当前值的显函数来计算下一个时间步的状态值。固定步长显式求解器以数学方式表示为：

其中x是状态；Dx是依赖于求解器的函数，用于估算状态导数；h是步长大小；n表示当前时间步。

从固定步长显式求解器的数学表达方式中不难看出欧拉公式，将如图1所示的输入u=1带入式（1.9）可以得到积分器模块的数学表达为：

联立式（1.10）与（1.11）可以得到积分器的输出为：

结合初始状态为0，这样就能解释如图二所示的积分器的输出值。

▲图二 积分器计算结果