MATLAB矩阵常微分方程数值解与数值仿真（连续时间多智能体系统）

在多智能体系统（MASs）论文数值仿真时，时常碰到连续系统的情况，其中参数也多为系数矩阵。本文主要对简单的连续系统微分方程（含系数矩阵）的MATLAB数值解求法进行归纳，并绘制仿真图像。同时也对一阶和高阶微分方程求法进行总结。

1. 可求解析解的微分方程

系统方程为一阶积分器的连续系统 $$\dot{x}=u,x\in R^n,$$
其中为了实现平均一致性，设计控制器 $u=-Lx$, $L\in R^{n\times n}$ 为对应图的拉普拉斯矩阵。
为了进行系统仿真,设 $n=4,L$ 如下，此处连续系统方程为
$$\dot{x}==-\left(\begin{array}{cccc}3& -1& -1& -1\\ -1& 2& -1& 0\\ -1& -1& 3& -1\\ -1& 0& -1& 2\end{array}\right)x$$
假设初值 $x_0=(10;0;5;10)$， 对其进行数值求解与仿真。

可见对于如上 $\dot{x}=-Lx$ 形式的微分方程，可用常微分方程的知识，先求出其解析解，再进行数值仿真。解析解为：
$$x=e^{-Lt}{x}_0.$$

MATLAB代码如下：

%%初始化设置
step=200;          %定义迭代步数
y=zeros(step,4);   %初始化矩阵来存储迭代值，用于作图
x0=[10; 0; 5; 10]; %微分方程初值
L=[3 -1 -1 -1;     %系数矩阵L
  -1  2 -1  0;
  -1 -1  3 -1;
  -1  0 -1  2];
for i=1:4         %迭代初值
    y(1,i)=x(i);
end

%%系统迭代步
for k=2:step      %迭代过程
    x=expm(-L*k/50)*x0;   %注意指数矩阵求解用函数expm()
    for i=1:4
        y(k,i)=x(i);
    end
    k=k+1;
end

%%仿真作图
t=1:step;
plot(t/50,y(t,1),t/50,y(t,2),t/50,y(t,3),t/50,y(t,4))
xlabel('t');
ylabel('x');
legend('x1','x2','x3','x4')

仿真图像如下：

2. MATLAB函数直接调用

在很多情况下，微分方程的解析解是很难求得的，所以才需要求其数值解并进行仿真。这里主要运用MATLAB自带函数ODE45（龙格库塔方法）来求解。
对于题设 1 中问题，我们将其转化为微分方程组，如下：
$$\{\begin{array}{c}\dot{x_1}=-[3x_1-x_2-x_3-x_4]\\ \dot{x_2}=-[-x_1+2x_2-x_3+0]\\ \dot{x_3}=-[-x_1-x_2+3x_3-x_4]\\ \dot{x_4}=-[-x_1+0-x_3+2x_4]\end{array}$$

定义m文件 $func.m$，将函数信息存入此m文件。

function dx=func(t,x)

%% 初始化有4个分量的dx
dx=zeros(4,1);

%% 微分方程
dx(1)=-(3*x(1)-x(2)-x(3)-x(4));    %dx(1)为x第一个分量的导数，下同
dx(2)=-(-x(1)+2*x(2)-x(3)+0);      %x(1)指x的第一个分量，下同
dx(3)=-(-x(1)-x(2)+3*x(3)-x(4));
dx(4)=-(-x(1)+0-x(3)+2*x(4));

如下代码对函数进行调用，并求数值解与仿真。

%% 参数初始化
startt=0;endd=10;      %区间开始和结束
x1=10;x2=0;x3=5;x4=10; %变量初值

%% ode45方法求解数值解
[t,x]=ode45(@func,[startt endd],[x1;x2;x3;x4]);

%% 仿真图像
plot(t,x(:,1),t,x(:,2),t,x(:,3),t,x(:,4))
xlabel('t');
ylabel('x');
legend('x1','x2','x3','x4')

仿真图像绘制如下：

3. 其他情形

3.1 一阶微分方程

简单的一阶微分方程，无需将微分方程写入m文件，只需在命令行的ODE45函数中加入微分方程即可。例如求解如下微分方程数值解
$$\dot{x}=x+t,$$

MATLAB代码如下：

%% 定义x初值
x0=0;            
%% ode5求微分方程数值解，其中[0 6]为求解区间
% @(t,x) x+t 为要求解的微分方程表示
[t,x]=ode45(@(t,x) x+t,[0 6],x0);
%% 绘制图像
plot(t,x);       

仿真图像显示如下：

3.2 高阶微分方程

高阶微分方程求数值解，基本思路是将其化为一阶微分方程组进行求解。例如对于如下二阶微分方程：
$$\begin{gathered} \ddot{x}+x^3+{\dot{x}}^3=0, \\ x(0)=1,\dot{x}(0)=0. \end{gathered}$$
我们将其化为两个一阶微分方程，如下：
$$\begin{gathered} {\dot{x}}_1=x_2, \\ {\dot{x}}_2=-{x_1}^3-{x_2}^3. \end{gathered}$$
其实，$x_1$表示原系统$x$, $x_2$表示原系统 $\dot{x}$. 化为一阶方程组后，便可采用 2 中介绍的方法，先定义M文件 $funcc.m$ 存储如上微分方程信息。MATLAB代码如下：

function dx=funcc(t,x)
dx=[x(2);
    -x(1)^3-x(2)^3];

命令行输入如下代码，调用m文件，并求解与绘图。

%% 给定x1和x2初值
x0=[1;0];

%% 使用ode45求解微分方程
[t,x]=ode45(@funcc,[0 20],x0);

%% 仿真作图（分别绘制了x和x的导函数图像）
plot(t,x(:,1),t,x(:,2))
xlabel('t');
ylabel('x');
legend('x_1','x_2')

仿真图像如下所示：

4. 总结

微分方程数值解求法总结：

1. 直接求出其解析解，便可计算其数值解，如情况1。
2. 不可求解析解的形式：
   2.1 一阶微分方程，如情况3.1。
   2.2 一阶微分方程组，如情况2。
   2.3 高阶微分方程，如情况3.2。

参考资料

[1] https://blog.csdn.net/qq_18820125/article/details/104727013
[2] https://www.zhihu.com/question/22378594