| MATLAB矩阵常微分方程数值解与数值仿真（连续时间多智能体系统） | 您所在的位置：网站首页 › matlab求含参积分 › MATLAB矩阵常微分方程数值解与数值仿真（连续时间多智能体系统） |
MATLAB矩阵常微分方程数值解与数值仿真（连续时间多智能体系统）
|
在多智能体系统（MASs）论文数值仿真时，时常碰到连续系统的情况，其中参数也多为系数矩阵。本文主要对简单的连续系统微分方程（含系数矩阵）的MATLAB数值解求法进行归纳，并绘制仿真图像。同时也对一阶和高阶微分方程求法进行总结。 目录 1. 可求解析解的微分方程2. MATLAB函数直接调用3. 其他情形3.1 一阶微分方程3.2 高阶微分方程 4. 总结参考资料 1. 可求解析解的微分方程系统方程为一阶积分器的连续系统 x ˙ = u , x ∈ R n , \dot{x}=u,~~x \in R^n, x˙=u, x∈Rn, 其中为了实现平均一致性，设计控制器 u = − L x u=-Lx u=−Lx, L ∈ R n × n L \in R^{n \times n} L∈Rn×n 为对应图的拉普拉斯矩阵。 为了进行系统仿真,设 n = 4 , L n=4,L n=4,L 如下，此处连续系统方程为 x ˙ = = − ( 3 − 1 − 1 − 1 − 1 2 − 1 0 − 1 − 1 3 − 1 − 1 0 − 1 2 ) x \dot{x}==-\begin{pmatrix} 3 & -1 & -1 & -1 \\ -1 & 2 & -1 & 0 \\ -1 & -1 & 3 & -1 \\ -1 & 0 & -1 &2 \\ \end{pmatrix}x x˙==−⎝⎜⎜⎛3−1−1−1−12−10−1−13−1−10−12⎠⎟⎟⎞x 假设初值 x 0 = ( 10 ; 0 ; 5 ; 10 ) x_0=(10;0;5;10) x0=(10;0;5;10)， 对其进行数值求解与仿真。 可见对于如上 x ˙ = − L x \dot{x}=-Lx x˙=−Lx 形式的微分方程，可用常微分方程的知识，先求出其解析解，再进行数值仿真。解析解为： x = e − L t x 0 . x=e^{-Lt}x_0. x=e−Ltx0. 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')仿真图像如下： 在很多情况下，微分方程的解析解是很难求得的，所以才需要求其数值解并进行仿真。这里主要运用MATLAB自带函数ODE45（龙格库塔方法）来求解。 对于题设 1 中问题，我们将其转化为微分方程组，如下： { x 1 ˙ = − [ 3 x 1 − x 2 − x 3 − x 4 ] x 2 ˙ = − [ − x 1 + 2 x 2 − x 3 + 0 ] x 3 ˙ = − [ − x 1 − x 2 + 3 x 3 − x 4 ] x 4 ˙ = − [ − x 1 + 0 − x 3 + 2 x 4 ] \left\{ \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} \right. ⎩⎪⎪⎨⎪⎪⎧x1˙=−[ 3x1−x2−x3−x4 ]x2˙=−[ −x1+2x2−x3+0 ]x3˙=−[ −x1−x2+3x3−x4 ]x4˙=−[ −x1+0−x3+2x4 ] 定义m文件 f u n c . m func.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')仿真图像绘制如下： 简单的一阶微分方程，无需将微分方程写入m文件，只需在命令行的ODE45函数中加入微分方程即可。例如求解如下微分方程数值解 x ˙ = x + t , \dot{x}=x+t, 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);仿真图像显示如下： 高阶微分方程求数值解，基本思路是将其化为一阶微分方程组进行求解。例如对于如下二阶微分方程： x ¨ + x 3 + x ˙ 3 = 0 , x ( 0 ) = 1 , x ˙ ( 0 ) = 0. \ddot{x}+x^3+\dot{x}^3=0,\\ x(0)=1,~\dot{x}(0)=0. x¨+x3+x˙3=0,x(0)=1, x˙(0)=0. 我们将其化为两个一阶微分方程，如下： x ˙ 1 = x 2 , x ˙ 2 = − x 1 3 − x 2 3 . \dot{x}_1=x_2,\\ \dot{x}_2=-{x_1}^3-{x_2}^3. x˙1=x2,x˙2=−x13−x23. 其实， x 1 x_1 x1表示原系统 x x x, x 2 x_2 x2表示原系统 x ˙ \dot{x} x˙. 化为一阶方程组后，便可采用 2 中介绍的方法，先定义M文件 f u n c c . m funcc.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')仿真图像如下所示： 微分方程数值解求法总结： 直接求出其解析解，便可计算其数值解，如情况1。不可求解析解的形式： 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 |
【本文地址】