[MATLAB作业二]求解线性代数问题

文章目录

0.0 题目
0.1 线性代数问题
- （1）用数值法和解析法求解下面的Sylvester方程，并验证得到的结果
- （2）试配平下面的化学反应方程式（不需要手算）

0.0 题目

0.1 线性代数问题

（1）用数值法和解析法求解下面的Sylvester方程，并验证得到的结果

[3−6−405142−24−63−673−13100−11004034]X+X[3−21−2−92−2−19]=[−21−14125−616−4−4−66−3]\begin{aligned} \begin{bmatrix} 3 & -6 & -4 & 0 & 5 \\ 1 & 4 & 2 & -2 & 4 \\ -6 & 3 & -6 & 7 & 3 \\ -13 & 10 & 0 & -11 & 0 \\ 0 & 4 & 0 & 3 & 4 \end{bmatrix} \boldsymbol{X}+\boldsymbol{X}\left[\begin{array}{rrr} 3 & -2 & 1 \\ -2 & -9 & 2 \\ -2 & -1 & 9 \end{array}\right]=\begin{bmatrix} -2 & 1 & -1 \\ 4 & 1 & 2 \\ 5 & -6 & 1 \\ 6 & -4 & -4 \\ -6 & 6 & -3 \end{bmatrix} \end{aligned} ⎣⎢⎢⎢⎢⎡31−6−130−643104−42−6000−27−11354304⎦⎥⎥⎥⎥⎤X+X⎣⎡3−2−2−2−9−1129⎦⎤=⎣⎢⎢⎢⎢⎡−2456−611−6−46−121−4−3⎦⎥⎥⎥⎥⎤

A=[3 -6 -4 0 5;1 4 2 -2 4;-6 3 -6 7 3;-13 10 0 -11 0;0 4 0 3 4]
B=[3 -2 1;-2 -9 2;-2 -1 9]
X0=[-2 1 -1;4 1 2;5 -6 1;6 -4 -4;-6 6 -3]

数值解

% 数值解
X = sylvester(A,B,-X0)

实验结果为：

X =-4.0569  -14.5128    1.56530.0356   25.0743   -2.74089.4886   25.9323   -4.41772.6969   21.6450   -2.88517.7229   31.9100   -3.7634

解析解

主要思想是设置关于ttt的指数形式，最后的求解包含ttt

A=sym(A)
B=sym(B)
X0=sym(X0)
syms t
% 解析解 dsolve方法
X(t)=any_matrix([5,3],'x',t)
x=dsolve(diff(X)==A*X+X*B,X(0)==X0)
X=[x.x11,x.x12,x.x13;x.x21,x.x22,x.x23,x.x31;...x.x32,x.x33,x.x41,x.x42,x.x43;x.x51,x.x52,x.x53]

% 解析解
X=simplify(expm(A*t)*X0*expm(B*t))
simplify(diff(X)-A*X-X*B)
subs(X,t,0)-X0

% 解析解 Kron方法(直积)
A0=kron(eye(3),A)+kron(B',eye(5));
X=expm(A0*t); x1=reshape(x,5,3)

其中，any_matrix函数如下，可以产生关于如x11(t),x12(t)x11(t),x12(t)x11(t),x12(t)这样的编号变量。

如

X(t) =[ x11(t), x12(t), x13(t)]
[ x21(t), x22(t), x23(t)]
[ x31(t), x32(t), x33(t)]
[ x41(t), x42(t), x43(t)]
[ x51(t), x52(t), x53(t)]

function A=any_matrix(nn,sA,varargin) % 生成任意矩阵
v=varargin; n=nn(1);
if length(nn)==1, m=n; else, m=nn(2); end
s=''; k=length(v); K=0; if n==1 || m==1, K=1; end
if k>0, s='('; for i=1:k, s=[s ',' char(v{i})]; ends(2)=[]; s=[s ')'];
end
for i=1:n, for j=1:m % 用循环结构逐个元素单独处理if K==0, str=[sA int2str(i),int2str(j) s];else, str=[sA int2str(i*j) s]; endeval(['syms ' str]); eval(['A(i,j)=' str ';']); % 指定相应的矩阵元素
end, end

遗憾的是，在处理四阶问题电脑可以顺利跑完，生成关于ete^tet的矩阵解，而对于五阶矩阵，电脑内存占用太多，上述方法运行失效，采用以下方法

function X=lyap3(A,B,C)
[nr,nc]=size(C);
A0=kron(A,eye(nc))+kron(eye(nr),B');
tryC1=C';X0=-inv(A0)*C1(:);X=reshape(X0,nc,nr)';
catcherror('矩阵奇异');
end

运行后结果如下,验证后为解析解

X = lyap3(sym(A),B,X0)ans =[ -434641749950/107136516451, -4664546747350/321409549353,   503105815912/321409549353]
[    3809507498/107136516451,  8059112319373/321409549353,  -880921527508/321409549353]
[ 1016580400173/107136516451,  8334897743767/321409549353, -1419901706449/321409549353]
[  288938859984/107136516451,  6956912657222/321409549353,  -927293592476/321409549353]
[  827401644798/107136516451, 10256166034813/321409549353, -1209595497577/321409549353]A*X+X*B+X0ans =[ 0, 0, 0]
[ 0, 0, 0]
[ 0, 0, 0]
[ 0, 0, 0]
[ 0, 0, 0]

（2）试配平下面的化学反应方程式（不需要手算）

KClO3+HCl→KCl+Cl2+H2O\begin{aligned} \mathrm{KClO}_{3}+\mathrm{HCl} \rightarrow \mathrm{KCl}+\mathrm{Cl}_{2}+\mathrm{H}_{2} \mathrm{O} \end{aligned} KClO3+HCl→KCl+Cl2+H2O

K4Fe(CN)6+KMnO4+H2SO4→KHSO4+Fe2(SO4)3+MnSO4+HNO3+CO2+H2\begin{aligned} \mathrm{K}_{4} \mathrm{Fe}(\mathrm{CN})_{6}+\mathrm{KMnO}_{4}+\mathrm{H}_{2} \mathrm{SO}_{4} \rightarrow \mathrm{KHSO}_{4}+\mathrm{Fe}_{2}\left(\mathrm{SO}_{4}\right)_{3}+\mathrm{MnSO}_{4}+\mathrm{HNO}_{3}+\mathrm{CO}_{2}+\mathrm{H}_{2} \end{aligned} K4Fe(CN)6+KMnO4+H2SO4→KHSO4+Fe2(SO4)3+MnSO4+HNO3+CO2+H2

配平方程式采用原子守恒方法，通过列等式来求得系数，目标函数为变量的系数和。

% KClO3 + HCl → KCl + Cl2 + H2O% KClO3 + 6HCl → 1KCl + 3Cl2 + 3H2O
% 变量包括
%     K C O H 以该顺序定义矩阵那么分别为
X1 = [1 1 3 0]
X2 = [0 1 0 1]
Y1 = [1 1 0 0]
Y2 = [0 2 0 0]
Y3 = [0 0 1 2]
A = [X1',X2']
b = [-Y1',-Y2',-Y3']
AA = [A,b]
%因此需要解的方程可以表示为
% AX1+BX2-CY1-DY2-EY3=0
% (A,B,-C.-D,-E)(X1,X2,X3,Y1,Y2,Y3)=0
% 因此在求解时候需要满足[x,fval] = linprog(ones(5,1),[],[],AA,zeros(4,1),ones(5,1))

实验结果如下，验证后该结果正确

Optimal solution found.
x =1.0000 6.0000 1.0000 3.0000 3.0000

在进行第二个方程配平时候，如果使用linprog会出现非整数解的情况，因此改为intlinprog，并将9个参数为整数的信息传给数组intcon

% 4K4Fe(CN)6 + 30KMnO4 + 82H2SO4 → 46KHSO4 + 2Fe2(SO4)3 + 30MnSO4 +24HNO3 + 24CO2 + 47H2% K4Fe(CN)6+KMnO4+H2SO4 → KHSO4+Fe2(SO4)3+MnSO4+HNO3+CO2+H2
% 变量包括
%     K F C N M H S O
X1 = [4 1 6 6 0 0 0 0]
X2 = [1 0 0 0 1 0 0 4]
X3 = [0 0 0 0 0 2 1 4]
Y1 = [1 0 0 0 0 1 1 4]
Y2 = [0 2 0 0 0 0 3 12]
Y3 = [0 0 0 0 1 0 1 4]
Y4 = [0 0 0 1 0 1 0 3]
Y5 = [0 0 1 0 0 0 0 2]
Y6 = [0 0 0 0 0 2 0 0]
A = [X1',X2',X3']
b = [-Y1',-Y2',-Y3',-Y4',-Y5',-Y6']AA = [A,b]
intcon = [1:9] %整数规划
[x,fval] = intlinprog(ones(9,1),intcon,[],[],AA,zeros(8,1),ones(9,1))

实验结果如下，验证后该结果正确

Optimal solution found.
x =4.0000 30.0000 82.0000 46.0000 2.0000 30.0000 24.0000 24.0000 47.0000