利用MATLAB求解一阶线性常系数非齐次微分方程组

主要步骤

本来想用在矩阵论期中开卷考试验证计算结果的，结果一个解方程组的题也没考…在一些学习网站白嫖惯了，也该分享一下知识，希望这里的code能对搜索到的人有帮助。
主要应用了matlab的符号矩阵，讲解不侧重推导，只侧重编程实现，详情可以看矩阵论相关的参考书。

1.问题形式

一阶线性常系数齐次微分方程组形式为：
$\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}$
将其改写为矩阵方程：
$\boldsymbol{x^{\prime}}=\frac{d\boldsymbol{x}}{d t}=\boldsymbol{Ax}+\boldsymbol{b}(t)$
其一般解为：
$\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$
其中$\boldsymbol{x}_0$为变量的初始值。

2.求矩阵函数

矩阵函数$e^{t\boldsymbol{A}}$的求解采用若当标准型法。
代码如下：
定义函数，输出为矩阵函数的矩阵

function [EtA,P,J] = my_expm( 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)!$的形式，存入数组备用。

%定义矩阵多项式的导数
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=\left[\begin{array}{rrr} -2 & 1 & 0 \\ -4 & 2 & 0 \\ 1 & 0 & 1 \end{array}\right]$
输出：

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_0$求解解得第一项

%先求矩阵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);

输入：
$\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$
输出：

ans =
              1
              1
 exp(t)*(t - 1)

转载：https://blog.csdn.net/Zhang_Engineer/article/details/105568535