使用Matlab作矩阵指数函数e^At和矩阵A

参考链接：

https://ww2.mathworks.cn/help/matlab/ref/expm.html?searchHighlight=expm&s_tid=doc_srchtitle

使用到的函数为expm 矩阵指数

语法

Y = expm(X)

说明

Y = expm (X) 计算 X 的矩阵指数。（即矩阵指数部分）

虽然不按此种方式计算，但是如果 X 包含一组完整的特征向量 P 和对应特征值 D

则 [P,D] = eig(X)

且（当矩阵A特征值互异的情形）expm (X) = P*diag(exp(diag(D)))*(P^(-1))

附：eig特征值和特征向量

矩阵A的相似变换

[P,D] = eig(A) 返回特征值的对角矩阵 D 和可逆矩阵P，其列是对应的右特征向量，使得 A*P = P*D。

矩阵A的对角化，可逆矩阵P和对角矩阵D。

(P^(-1))*A*P=D

A*P = P*D

示例

求矩阵A=

   -2  0  0

0 -3  1

0  0 -3

的矩阵指数函数e^At

输入命令行

a=[-2 0 0;0 -3 1;0 0 -3];   %定义矩阵a

>> syms t;               %定义变量t

>> expm(a*t)             %利用expm函数计算转移矩阵

得到结果：

ans =

[ exp(-2*t),         0,           0]

[         0, exp(-3*t),   t*exp(-3*t)]

[         0,       0,     exp(-3*t)]

通过拉氏反变换来求e^At

输入命令行：

a=[-2 0 0;0 -3 1;0 0 -3];      %定义矩阵a

>> syms s t;                %定义vi按量s,t

>> G=inv(s*eye(size(a))-a);   %(sI-a)的逆

>> ilaplace(G)              %拉氏反变换

ans =

[ exp(-2*t),         0,           0]

[         0, exp(-3*t),   t*exp(-3*t)]

[         0,       0,     exp(-3*t)]

附：

inv逆矩阵的求取

Y = inv(X) 计算方阵X的逆矩阵。

X^(-1) 等效于inv(X)。

size数组的大小

sz = size(A) 返回一个行向量，其元素包含 A 的相应维度的长度。例如，如果 A 是一个 3×4 矩阵，则 size(A) 返回向量 [3 4]。

eye单位矩阵

I = eye(n,m) 返回一个主对角线元素为 1 且其他位置元素为 0 的 n×m 矩阵。

ilaplace拉氏反变换

F(t)=e^At

F’(t)=At*(e^At)      这是矩阵对标量的求导

令t=0, F’(0)=A

示例

已知矩阵F(t)=[  1   (1/2)*(1-exp(-2*t))   ;   0    exp(-2*t)   ]

(一个2x2矩阵)

求A矩阵

在Matla输入命令行：

syms t;

>> a=[1 (1/2)*(1-exp(-2*t));0 exp(-2*t)];

>> da=diff(a)

得到结果：

da =

[ 0,    exp(-2*t)]

[ 0, -2*exp(-2*t)]

附：

diff

差分和近似导数

Y = diff(X)

参考文档：

https://ww2.mathworks.cn/help/matlab/ref/diff.html?searchHighlight=diff&s_tid=doc_srchtitle

x‘=Ax+Bu

y=Cx+Du

A=[0 1;0 0],B=[0 ; 1],C=[1 0];D=0;

x'=[0 1;0 0]x+[0 ; 1]u

y=[1 0]x

初始状态x(0)=[1 ; 1]，输入u(t)是单位阶跃响应

输入命令行

A=[0 1;0 0];B=[0;1];

C=[1 0];D=0;            %定义系统矩阵

G=ss(A,B,C,D);          %建立状态空间描述的系统模型

x0=[1;1];               %初始状态

syms s t;               %定义变量

G0=inv(s*eye(size(A))-A); %(sI-a)的逆

x1=ilaplace(G0)*x0;     %零输入响应x1

G1=inv(s*eye(size(A))-A)*B;

x2=ilaplace(G1/s)    %1/s为积分

x=x1+x2;

y=C*x;

得出结果：

x1 =

 t + 1

     1

x2 =

 t^2/2

     t

x =

 t^2/2 + t + 1

         t + 1

y =

t^2/2 + t + 1