MATLAB解矩阵微分方程

1 2 3 4 5 6
[ t , y ] = r u n g e K u t t a 4 ( @odefunc , 0 , 1 0 , [ 0 0 ] , 0 . 1 ) ; % 调 用 4 阶 龙 格 库 塔 法 figure (1) p l o t ( t , y ( : , 1 ) ) %绘 制 质 量 块 位 移 随 时 间 变 化 的 图 形 xlabel ( ’ 时间 t / s ’ ) y l a b e l ( ’ 位 移 y / m’ ) t i t l e ( ’ 弹簧质量块儿阻尼器系统 ’) 程序运行的结果如图 2
1 2 3 4 5
f u n c t i o n [ dydt ] = odefunc ( t , y ) %ODEFUNC 弹 簧 质 量 块 阻 尼 器 系 统 动 力 学 微 分 方 程 % % % y (1)： 状态变量 —— 位移 y (2)： 状态变量 —— 速度 dydt ( 1 ) ： 状 态 变 量 为 y时 ， 第 一 个 状 态 变 量 y (1) 的变化率 2
MATLAB 求解矩阵微分方程
鲁鹏 北京理工大学宇航学院 2018 年 06 月 25 日
本文先简要介绍了矩阵的积分和微分的定义； 接着通过弹簧-质量块儿-阻尼器系统的例子， 阐述了求解常微分方程数值解的常用方法；最后，在前两部分的基础上，总结了用 MATLAB ˙ (t) = f (X (t), t) 矩阵微分方程数值解的方法。 求解形如 X
ft = linspace (0 ,30 ,300); f = 5* s i n ( 1 . 2 5 * s q r t ( k1 / m1) * f t ) ; % ode45 求 解 ， 也 可 以 使 用 4 阶 龙 格 库 塔 法 ， 但 是 要 注 意 参 数 格 式 不 同 。 [ t , x ] = ode45 (@( t , x ) o d e f u n c 1 ( t , x , A, B , f t , f ) , [ 0 3 0 ] , [ 0 0 0 0 0 0 ] ) ; figure (1) p l o t ( t , x ( : , 1 ) ) %绘 制 质 量 块 1 位 移 随 时 间 变 化 的 图 形 h o l d on p l o t ( t , x ( : , 2 ) ) %绘 制 质 量 块 2 位 移 随 时 间 变 化 的 图 形 p l o t ( t , x ( : , 3 ) ) %绘 制 质 量 块 3 位 移 随 时 间 变 化 的 图 形 xlabel ( ’ 时间 t / s ’ )
k3 m3 k3 −m 3
0 0 x1 0 x2 0 1 x3 0 + 1 0 x4 m1 0 x5 0 x6 0
0 0 0 0
1 m2
0

0 f1 0 f2 0 f 3 0
f u n c t i o n [ x , y ] = r u n g e K u t t a 4 ( u f u n c , a , b , y0 , h ) % 参 数 表 顺 序 依 次 是 微 分 方 程 组 的 函 数 句 柄 ufunc ， 时 间 起 点 a ， % 时 间 终 点 b ， 初 始 值 向 量 y0 ， 步 长 h （ 参 数 形 式 参 考 了 ode45 函 数 ） n = f l o o r ( ( b − a ) / h );% 求 步 数 x = ones ( n , 1 ) ; x ( 1 ) = a ;% 时 间 起 点 y ( : , 1 ) = y0;% 赋 初 值 ， 可 以 是 向 量 ， 但 是 要 注 意 维 数 for i i = 1: n x ( i i + 1) = x ( i i ) + h ; k1 = u f u n c ( x ( i i ) , y ( : , i i ) ) ; k2 = u f u n c ( x ( i i ) + h / 2 , y ( : , i i ) + h * k1 / 2 ) ; k3 = u f u n c ( x ( i i ) + h / 2 , y ( : , i i ) + h * k2 / 2 ) ; k4 = u f u n c ( x ( i i ) + h , y ( : , i i ) + h * k3 ) ; y ( : , i i + 1 ) = y ( : , i i ) + h * ( k1 + 2 * k2 + 2 * k3 + k4 ) / 6 ; end y = y ’; end 使用龙格库塔法求解方程 1时，可以先新建一个函数文件 odefunc.m，该函数文件包含了常微 分方程 1的所有信息。odefunc.m 的内容如下：
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
m1 = 1 ; m2 = 1 ; m3 = 1 ; k1 = 1 ; k2 = 1 ; k3 = 1 ; A = [0 0 0 −(k1+k2 ) / m1 k2 / m2 0 B = [0 0 0 1 / m1 0 0 0 0 0 0 1 / m2 0 0 0 0 k2 / m1 −(k2+k3 ) / m2 k3 / m3 0; 0; 0; 0; 0; 1 / m3 ] ; 0 0 0 0 k3 / m2 −k3 / m3 1 0 0; 0 1 0; 0 0 1; 0 0 0; 0 0 0; 0 0 0];
4 5 6 7 8
% y ’ = f ( t , y ) 从 t 0 =0 到 t f =10 的 积 分 ， [ 0 0 ] 为 初 始 条 件 y0 ) p l o t ( t , y ( : , 1 ) ) %绘 制 质 量 块 位 移 随 时 间 变 化 的 图 形 xlabel ( ’ 时间 t / s ’ ) y l a b e l ( ’ 位 移 y / m’ ) t i t l e ( ’ 弹簧质量块儿阻尼系统 ’) 使用 MATLAB 的常微分方程求解器 ode45 求解的结果如下图 3：
1 m3
(2)
0 0 0 4
0
图 4: 三自由度弹簧质量块儿系统 将系统状态方程写成 MATLAB 函数文件 odefunc1.m，odefunc1.m 的内容如下：
1 2 3 4 5 6
f u n c t i o n [ dx ] = o d e f u n c 1 ( t , x , A, B , f t , f ) %ODEFUNC1 三 自 由 度 弹 簧 质 量 块 儿 系 统 状 态 方 程 % 此处显示详细说明 f = interp1 ( ft , f , t ) ; % 对数据 ( ft , f ) 进行插值以获取 t 时刻的 f 值 dx = A*x+B* [ f f f ] ’ ; end 取 k1 = k2 = k3 = 1，m1 = m2 = m3 = 1，f1 = f2 = f3 = 5 sin(1.25t)，求解方程 2数值解的 MATLAB 程序如下：
0.35
0.3
0.25
0.2
0.15
0.1
0.05
0 1 2 3 4 5 6 7 8 9 10
图 3: ode45 求解的质量块儿位移随时间变化的图像
3
矩阵微分方程数值解
˙ (t) = AG(t) + G(t)AT + BR−1 B T 是常见的两种矩阵微分方程， x ˙ (t) = Ax(t) + Bu(t) 和 G
0.35
0.3
0.25
0.2
0.15
0.1
0.05
0 0 2 4 6 8 10
图 2: 4 阶龙格库塔法求解的质量块儿位移随时间变化的图像
2.2
常微分方程求解器 ode45
使用 MATLAB 的常微分方程求解器 ode45 也可以求解上一小节质量块儿位移随，MAT-
LAB 代码如下：
1 2 3
figure [ t , y ] = ode45 ( @odefunc , [ 0 1 0 ] , [ 0 0 ] ) ; % ode45 ( @odefunc 是 函 数 o d e f u n c ( ) 的 句 柄 ， [ 0 1 0 ] 表 示 求 微 分 方 程 组 3
6 7 8 9 10
%
dydt ( 2 ) ： 状 态 变 量 为 y时 ， 第 二 个 状 态 变 量 y (2) 的变化率
dydt = zeros ( 2 , 1 ) ; dydt (1) = y ( 2 ) ; d y d t ( 2 ) = −5 * y ( 1 ) − 0 . 7 * y ( 2 ) + 1 ; end 接着新建一个脚本文件， 文件名任意， 这里使用的文件名是 main.m， 脚本文件 main.m 的内容如 下代码框中所示。这里需要注意的是必须要将两个函 0 数文件 odefunc.m 和 rungejielonggekutamain.m 放在同一个文件夹中。
这里以这两类矩阵微分方程求解为例，总结 MATLAB 求解这矩阵微分方程数值解的方法。
3.1
第一个例子
图 4是一个三自由度弹簧质量块儿系统 k1 ,k2 ,k3 分别是三个弹簧的弹性系数,m1 ，m2 ，m3
分别是三个质量块儿的质量，f1 ,f2 ,f3 分别是给三个质量块儿施加的外力，外力是时间 t 的函 数。对于该三自由度弹簧质量块儿系统，取质量块儿的位移 y1 ，y2 ，y3 和速度 v1 ，v2 ，v3 为 状态变量，则系统的状态方程 [1] 可以表示为： 0 0 0 1 0 x ˙1 x 0 0 0 1 ˙ 0 2 x 0 0 0 0 ˙ 3 0 = k1 +k2 k 2 x 0 0 0 m1 ˙ 4 − m1 k2 k3 k2 +k3 ˙ 5 m2 − m2 0 0 x m2 x ˙6 0
t0 tΒιβλιοθήκη Baidu m ×n
即对矩阵积分，就是对矩阵每个元素积分。
2
常微分方程数值解
求常微分方程数值解常用的方法有 4 阶龙格库塔法，也可以使用 MATLAB 的常微分方程
求解器 ode45。
图 1: 弹簧-质量块儿-阻尼器系统结构图
1
此处以图 1所示弹簧-质量块-阻尼器系统为例。总结用 MATLAB 求解常微分方程数值解 的方法。质量块儿质量 m = 1kg ，弹簧弹性系数 k = 5，阻尼器阻尼系数 f = 0.7，输入的 力 F = 10N (t ⩾ 0)，则动力学方程为 y ¨(t) + 0.7y ˙ (t) + 5y (t) = 1。取位移 y1 = y (t) 和速度 y2 = y ˙ (t) 为状态变量，将此方程重写为一阶常微分方程 y ˙ 1 = y2 y ˙ 2 = −5y1 − 0.7y2 + 1
(1)
2.1
4 阶龙格库塔法
《数值分析》中有龙格库塔法的原理，龙格库塔法 MATLAB 程序如下：
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
% RUNGEKUTTA4 4 阶 龙 格 库 塔 法 % 作 者 ： h y o w i n n e r (MATALB 中 文 论 坛 用 户 ) % 创建日期 ：未 知 % 修 改 人：鲁 鹏( 北京理工大学宇航学院 ) % 修 改 日 期 ： 2 0 18 年 6 月 24 日 % 参 考 网 站 ： h t t p : / / www. i l o v e m a t l a b . cn / t h r e a d − 42808 − 1 − 1. h t m l % 版 本 ：2.0
1
矩阵的积分和微分
( )
m×n
Definition 1 若矩阵 A(t) = (aij (t))m×n 的每一个元素 aij (t) 是变量 t 的可微函数，则称 A(t) 可微，其导数定义为 dA = A′ (t) = dt daij dt
即对矩阵求导，就是对矩阵每个元素求导。由此出发，函数可以定义高阶导数。 Definition 2 若矩阵 A(t) = (aij (t))m×n 的每一个元素 aij (t) 都是区间 [t0 , t1 ] 上的可积函数， 则 称 A(t) 在区间 [t0 , t1 ] 上可积，并定义 A(t) 在区间 [t0 , t1 ] 上的积分为 (∫ t1 ) ∫ t1 A(t)dt = aij (t)dt