主要步骤

本来想用在矩阵论期中开卷考试验证计算结果的，结果一个解方程组的题也没考…在一些学习网站白嫖惯了，也该分享一下知识，希望这里的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:nif J(i-1,i) == 0J_{num} = J(pre:i-1,pre:i-1);pre = i;num = num + 1;if i == n J_{num} = J(pre,pre);endendif i == n && J(i-1,i) == 1J_{num} = J(pre:n,pre:n);end
end

搜索尺寸最大的若当块，因为特征多项式的导数次数由最大的若当块尺寸决定，这个尺寸可以和上一次遍历合并，这样一次遍历就可以得到所有的若当块和最大若当块尺寸。

%搜索最大维度
mSize = 0;
for i = 1:numif size(J_{i},1) > mSizemSize = 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-1f_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 >= 2for j = 1:n1-1%若当块维度索引etJ = etJ + diag(subs(f_J(j+1),x,J_{i}(1,1))*ones(1,n1-j),j);endendEtJ = 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 =11exp(t)*(t - 1)

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

[每日一氵]求解一阶线性常系数微分方程组
求解一阶线性常系数微分方程组关键字有这么多: 一阶线性常系数微分方程组直接给例子吧: d x 1 d t = x 2 d x 2 d t = x 3 d x 3 d t = − 6 x 1 ...
matlab求解全微分函数,利用MATLAB求解微分方程的方法探索
引言科学问题和工程问题经常需要求取微分方程的解,MATLAB 的强大数值运算和符号运算能力,能够方便地进行各种解析运算,是方便实用.功能强大的数学软件之一. 1线性微分方程求解 1.1线性常微分方程 ...
线性规划问题的数学建模matlab,数学建模讲座之三——利用Matlab求解线性规划问题（linprog函数）.ppt...
数学建模讲座之三--利用Matlab求解线性规划问题(linprog函数) 利用Matlab求解线性规划问题线性规划是一种优化方法,Matlab优化工具箱中有现成函数linprog对如下式描述的LP ...
利用matlab求解线性规划,数学建模讲座之三利用matlab求解线性规划问题（linprog函数）...
数学建模讲座之三利用matlab求解线性规划问题(linprog函数) 利用利用 Matlab求解线性规划问题求解线性规划问题河北科技河北科技大学大学*第第 1页页利用 Matlab求解线性规划问题 ...
运筹学matlab实验报告,运筹学上机实验报告利用Matlab求解整数线性规划
四川师范大学数学与软件科学学院运筹学上机实验报告. 学期:__2011_至__2012__ 第___一__ 学期 2011年11月9日课程名称:__ 运筹学 ________ 专业:_信息与计算 ...
python求解一阶线性偏微分方程通解举例
python求解一阶线性偏微分方程的通解举例 Python求解偏微分方程也是其一个应用方面,下面举例说明. 一.问题: 求一阶线性偏微分方程 x ∂ f ( x , y ) ∂ x − y ∂ f ( ...
在matlab中实现累乘,如何利用matlab设计一个线性相位FIR带通滤波器，并在FPGA上实现...
设计要求利用matlab设计一个线性相位FIR带通滤波器,并在FPGA上实现. 1.滤波器指标:过渡带带宽分别为100~300HZ,500~700HZ,阻带允许误差为0.02,通带允许误差为0.01 ...
如何利用matlab求解方程
如何利用matlab求解方程 1. 前言作为三大数学软件之一,matlab在数值计算方法的能力首屈一指.求解方程是工科学习和工程计算中最基础.最常见的问题.掌握利用现代化工具求解方程的方法对于 ...
matlab 线性时不变规律,MATLAB实验——运用MATLAB求解和线性时不变系统要点详解.docx...
MATLAB实验报告课程名称 MATLAB程序设计实验日期 2015 年 05 月 18 日学生姓名学号班级实验名称运用MATLAB求解和分析线性时不变系统实验仪器MATLAB7.1 Win ...

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

用矩阵函数求解一阶线性常系数非齐次微分方程组

主要步骤

1.问题形式

2.求矩阵函数

3.代入矩阵A的指数函数得最终解

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

最新文章

热门文章