用矩阵函数求解一阶线性常系数非齐次微分方程组
主要步骤
本来想用在矩阵论期中开卷考试验证计算结果的，结果一个解方程组的题也没考…在一些学习网站白嫖惯了，也该分享一下知识，希望这里的code能对搜索到的人有帮助。
主要应用了matlab的符号矩阵，讲解不侧重推导，只侧重编程实现，详情可以看矩阵论相关的参考书。
1.问题形式
一阶线性常系数齐次微分方程组形式为：
d ξ 1 d t = a 11 ξ 1 + a 12 ξ 2 + ⋯ + a 1 n ξ n + β 1 ( t ) d ξ 2 d t = a 21 ξ 1 + a 22 ξ 2 + ⋯ + a 2 n ξ n + β 2 ( t ) ⋮ d ξ n d t = a n 1 ξ + a n 2 ξ 2 + ⋯ + a m ξ n + β n ( t ) \begin{array}{rl} \frac{\mathrm{d} \xi_{1}}{\mathrm{d} t} & =a_{11} \xi_{1}+a_{12} \xi_{2}+\cdots+a_{1 n} \xi_{n} +\beta_1(t) \\ \\ \frac{\mathrm{d} \xi_{2}}{\mathrm{d} t} & =a_{21} \xi_{1}+a_{22} \xi_{2}+\cdots+a_{2 n} \xi_{n} +\beta_2(t) \\ & \vdots \\ \frac{\mathrm{d} \xi_{n}}{\mathrm{d} t} & =a_{n 1} \xi+a_{n 2} \xi_{2}+\cdots+a_{m} \xi_{n}+\beta_n(t) \end{array} dtdξ1dtdξ2dtdξn=a11ξ1+a12ξ2+⋯+a1nξn+β1(t)=a21ξ1+a22ξ2+⋯+a2nξn+β2(t)⋮=an1ξ+an2ξ2+⋯+amξn+βn(t)
将其改写为矩阵方程：
x ′ = d x d t = A x + b ( t ) \boldsymbol{x^{\prime}}=\frac{d\boldsymbol{x}}{d t}=\boldsymbol{Ax}+\boldsymbol{b}(t) x′=dtdx=Ax+b(t)
其一般解为：
x ( t ) = e ( t − t 0 ) A x 0 + e t A ∫ t 0 t e − s A b ( s ) d s \boldsymbol{x}(t)=\mathbf{e}^{\left(t-t_{0}\right) \boldsymbol{A}} \boldsymbol{x}_{0}+\mathrm{e}^{t\boldsymbol{A}} \int_{t_{0}}^{t} \mathrm{e}^{-s\boldsymbol{A}} \boldsymbol{b}(s) \mathrm{d} s x(t)=e(t−t0)Ax0+etA∫t0te−sAb(s)ds
其中 x 0 \boldsymbol{x}_0 x0为变量的初始值。
2.求矩阵函数
矩阵函数 e t A e^{t\boldsymbol{A}} etA的求解采用若当标准型法。
代码如下：
定义函数，输出为矩阵函数的矩阵
function [EtA,P,J] = my_expm( A )
首先求若当标准型，这里注意一定要将矩阵 A A A转化为符号矩阵形式，标准型和特征向量会以分数和根号形式表示，并将所有若当块用元组（cell）存储
syms t x%x表示奇异值
[P,J] = jordan(sym(A));%先求若当标准型
n = size(A,1);
J_ = { };
%num存储若当块的个数
num = 1;pre = 1;
for i = 2:n
if J(i-1,i) == 0
J_{ num} = J(pre:i-1,pre:i-1);
pre = i;
num = num + 1;
if i == n
J_{ num} = J(pre,pre);
end
end
if i == n && J(i-1,i) == 1
J_{ num} = J(pre:n,pre:n);
end
end
搜索尺寸最大的若当块，因为特征多项式的导数次数由最大的若当块尺寸决定，这个尺寸可以和上一次遍历合并，这样一次遍历就可以得到所有的若当块和最大若当块尺寸。
%搜索最大维度
mSize = 0;
for i = 1:num
if size(J_{ i},1) > mSize
mSize = size(J_{ i},1);
end
end
在这里定义特征多项式和其 n − 1 n-1 n−1次导数除以 ( n − 1 ) ! (n-1)! (n−1)!的形式，存入数组备用。
%定义矩阵多项式的导数
f_J = [exp(t*x)];
for j = 1:mSize-1
f_J = [f_J,diff(exp(t*x),x,j)/prod(1:j)];
end
f_J = simplify(f_J);
对每一个若当块求特征函数矩阵，将所有特征函数矩阵对角合并到一起。
EtJ = [];
%blkdiag对角合并
for i = 1:num%若当块索引
n1 = size(J_{ i},1);
etJ = subs(f_J(1),x,J_{ i}(1,1))*diag(ones(1,n1));
if n1 >= 2
for j = 1:n1-1%若当块维度索引
etJ = etJ + diag(subs(f_J(j+1),x,J_{ i}(1,1))*ones(1,n1-j),j);
end
end
EtJ = blkdiag(EtJ,etJ);%若当块合并
end
EtJ = simplify(EtJ);
EtA = P*EtJ*(inv(P));
输入：
A = [ − 2 1 0 − 4 2 0 1 0 1 ] A=\left[\begin{array}{rrr} -2 & 1 & 0 \\ -4 & 2 & 0 \\ 1 & 0 & 1 \end{array}\right] A=⎣⎡−2−41120001⎦⎤
输出：
ans =
[ 1 - 2*t, t, 0]
[ -4*t, 2*t + 1, 0]
[ 2*t - exp(t) + 1, exp(t) - t - 1, exp(t)]
3.代入矩阵A的指数函数得最终解
定义函数
function [x,EtA] = my_solve_equation( A,b,co,to )
%输入：A系数矩阵，b非齐次项（齐次为0）,co初始值,to初始时刻
%一阶线性常系数非齐次\齐次微分方程组
将 t t t替换为 t − t 0 t-t_0 t−t0求解解得第一项
%先求矩阵A指数函数
[EtA,~,~] = my_expm( A,1 );
syms t s
Firs = subs(EtA,t,t-to)*co;
求解第二项矩阵积分项
b = subs(b,t,s);%一步符号替换
E_tA = subs(EtA,t,-s);
Secn = EtA*int(E_tA*b,s,to,t);
x = simplify(Firs + Secn);
输入：
A = [ − 2 1 0 − 4 2 0 1 0 1 ] b ( t ) = [ 1 2 e t − 1 ] x ( 0 ) = [ 1 1 − 1 ] t 0 = 0 \boldsymbol{A}=\left[\begin{array}{rrr} -2 & 1 & 0 \\ -4 & 2 & 0 \\ 1 & 0 & 1 \end{array}\right]\quad \boldsymbol{b}(t)=\left[\begin{array}{rrr} 1\\ 2\\ e^{t}-1 \end{array}\right]\quad \boldsymbol{x}(0)=\left[\begin{array}{rrr} 1\\ 1\\ -1 \end{array}\right]\quad t_0 = 0 A=⎣⎡−2−41120001⎦⎤b(t)=⎣⎡12et−1⎦⎤x(0)=⎣⎡11−1⎦⎤t0=0
输出：
ans =
1
1
exp(t)*(t - 1)
