第1章矩阵与行列式【矩阵与行列式简介】在计算机日益发展的今天，线性代数起着越来越重要的作用。线性代数起源于解线性方程组的问题，而利用矩阵来求解线性方程组的Gauss消元法至今仍是十分有效的计算机求解线性方程组的方法。矩阵是数学研究和应用的一个重要工具，利用矩阵的运算及初等变换可以解决求解线性方程组等问题。特殊的矩阵方阵的数字特征之一是方阵的行列式，使用行列式可以描述方阵的一些重要的性质。通过计算行列式可求逆矩阵，n个

第1章矩阵与行列式未知量n个方程的线性方程组的惟一解等问题。向量也是研究矩阵的有力工具，可通过向量组的秩来定义矩阵的秩。向量与矩阵、行列式都是线性代数的重要基本概念，它们是建立线性方程组的解的构造理论与系统求解方法的三个基本工具。

第1章矩阵与行列式验证性实验实验一矩阵的运算【实验目的】 1．理解矩阵、逆矩阵的概念
第1章矩阵与行列式验证性实验实验一矩阵的运算【实验目的】 1．理解矩阵、逆矩阵的概念 2．掌握矩阵的线性运算、乘法、转置、逆、方阵的幂的运算【实验要求】理解矩阵赋值命令、符号变量说明syms、加法+、乘法*、转置’、逆矩阵inv、方阵的幂^等命令

第1章矩阵与行列式【实验内容】 1．已知下列矩阵：（1），；（2），．计算，，，，，，．

第1章矩阵与行列式【实验过程】 1．（1）>> A=[3 1 1;2 1 2;1 2 3]; >> B=[1 1 -1;2 -1 0;1 0 1]; >> C=A+B 运行结果： C =

第1章矩阵与行列式 >> AB=A*B 运行结果： AB = >> D=6*A D =

第1章矩阵与行列式 >> sym c; >> cA=c*A 运行结果： cA = [ 3*c, c, c] [ 2*c, c, 2*c] [ c, 2*c, 3*c] >> F=A' F =

第1章矩阵与行列式 >> G=inv(A) 运行结果： G = 1/4 1/4 -1/ /4 5/4 -1/4 >> H=A^5 H =

第1章矩阵与行列式（2）>> A=sym('[a b;c d]');
第1章矩阵与行列式（2）>> A=sym('[a b;c d]'); >> B=sym('[1 a;1 b]'); >> C=A+B 运行结果： C = [ a+1, b+a] [ c+1, d+b] >> AB=A*B AB = [ b+a, a^2+b^2] [ c+d, c*a+d*b]

第1章矩阵与行列式 >> D=6*A 运行结果： D = [ 6*a, 6*b] [ 6*c, 6*d] >> sym c; >> cA=c*A cA = [ c*a, c*b] [ c^2, c*d]

第1章矩阵与行列式 >> F=A' 运行结果： F = [ conj(a), conj(c)] [ conj(b), conj(d)] % conj为复数共轭即 >> G=inv(A) G = [ d/(a*d-c*b), -b/(a*d-c*b)] [ -c/(a*d-c*b), a/(a*d-c*b)] 即．

第1章矩阵与行列式实验二矩阵的初等变换【实验目的】 1．理解矩阵初等变换的概念 2．掌握矩阵的初等变换及用初等变换求矩阵的逆矩阵
第1章矩阵与行列式实验二矩阵的初等变换【实验目的】 1．理解矩阵初等变换的概念 2．掌握矩阵的初等变换及用初等变换求矩阵的逆矩阵【实验要求】掌握矩阵的表示、符号变量说明syms、逆矩阵inv等命令

第1章矩阵与行列式【实验内容】 1．已知矩阵，求对矩阵实施如下的初等变换后所得矩阵。矩阵的第2行乘以m；
第1章矩阵与行列式【实验内容】 1．已知矩阵，求对矩阵实施如下的初等变换后所得矩阵。矩阵的第2行乘以m；矩阵的第3列的n倍加到第1列上去；矩阵的第1行与第2行交换。 2．已知矩阵，提取矩阵的第2、3、4行与第3、4列的元素构成矩阵B，提取矩阵的第2、3、4行与第1、4列的元素构成矩阵C．

第1章矩阵与行列式 3. 用初等变换求矩阵的逆矩阵。 4．已知，，且，求．

第1章矩阵与行列式【实验过程】 1．（1）>> syms m;
第1章矩阵与行列式【实验过程】 1．（1）>> syms m; >>A=sym('[a b c d;e f g h;i j k l]'); >>A(2,:)=m*A(2,:) 运行结果： A = [ a, b, c, d] [ m*e, m*f, m*g, m*h] [ i, j, k, l] （2）>> syms n; >>A(:,1)=A(:,1)+n*A(:,3) [ a+n*c, b, c, d] [ e+n*g, f, g, h] [ i+n*k, j, k, l]

第1章矩阵与行列式（3）>> A=sym('[a b c d;e f g h;i j k l]');
第1章矩阵与行列式（3）>> A=sym('[a b c d;e f g h;i j k l]'); >>A([2,1],:)=A([1,2],:) 运行结果： A = [ e, f, g, h] [ a, b, c, d] [ i, j, k, l] 2．>> A=[ ; ; ; ]; >>B=A(2:4,3:4) B =

第1章矩阵与行列式 >>C=A(2:end,[1,4]) 运行结果： C = 5 8 9 12 13 16
第1章矩阵与行列式 >>C=A(2:end,[1,4]) 运行结果： C = 3．>> A=[0 1 2;1 1 4;2 -1 0]; >>E=eye(3); >>B=[A,E] B =

第1章矩阵与行列式 >> B([1 2],:)=B([2 1],:) 运行结果： B = 1 1 4 0 1 0
第1章矩阵与行列式 >> B([1 2],:)=B([2 1],:) 运行结果： B = >> B(3,:)=B(3,:)-2*B(1,:)

第1章矩阵与行列式 >> B(3,:)=B(3,:)+3*B(2,:) 运行结果： B = 1 1 4 0 1 0
第1章矩阵与行列式 >> B(3,:)=B(3,:)+3*B(2,:) 运行结果： B = >> B(2,:)=B(2,:)+B(3,:); >>B(1,:)=B(1,:)+2*B(3,:)

第1章矩阵与行列式 >> B(1,:)=B(1,:)-B(2,:); >>B(3,:)=-1/2*B(3,:)
第1章矩阵与行列式 >> B(1,:)=B(1,:)-B(2,:); >>B(3,:)=-1/2*B(3,:) 运行结果： B = / /2 4．>> A=[1 0 1;-1 1 1;2 -1 1]; >>B=[1 1; 0 1;-1 0]; >>X=inv(A)*B X =

第1章矩阵与行列式实验三 Gauss消元法【实验目的】掌握解线性方程组的Gauss消元法
第1章矩阵与行列式实验三 Gauss消元法【实验目的】掌握解线性方程组的Gauss消元法【实验要求】掌握矩阵赋值命令、初等变换相关命令、简化矩阵为阶梯形式rref等命令【实验内容】 1．用Gauss消元法解线性方程组：（1）；（2）．

第1章矩阵与行列式【实验过程】 1．（1）解法一：Gauss消元法．
第1章矩阵与行列式【实验过程】 1．（1）解法一：Gauss消元法． >>A=[ ; ; ; ] ; >>A(2,:)=A(2,:)-A(1,:); >>A(3,:)=A(3,:)-2*A(1,:); >>A(4,:)=A(4,:)-A(1,:) 运行结果： A = >> A([2,3],:)=A([3,2],:)

第1章矩阵与行列式 >>A(2,:)=(-1)*A(2,:); >>A(3,:)=1/2*A(3,:) 运行结果：
第1章矩阵与行列式 >>A(2,:)=(-1)*A(2,:); >>A(3,:)=1/2*A(3,:) 运行结果： A = >> A(4,:)=A(4,:)-A(3,:); >>A(1,:)=A(1,:)-A(3,:); >>A(2,:)=A(2,:)-A(3,:)

第1章矩阵与行列式 >>A(1,:)=A(1,:)-2*A(2,:) 运行结果： A = 1 0 0 3 0 1 0 2
第1章矩阵与行列式 >>A(1,:)=A(1,:)-2*A(2,:) 运行结果： A = 由上可知，方程组有惟一解．解法二：>> A=[ ; ; ; ]; >>A=rref(A) 由上可知，结果同解法一。

第1章矩阵与行列式（2）解法一：Gauss消元法．
第1章矩阵与行列式（2）解法一：Gauss消元法． >> A=[ ; ; ] ; >> A([1,3],:)=A([3,1],:) 运行结果： A = >> A(2,:)=A(2,:)+A(1,:); >> A(3,:)=A(3,:)-2*A(1,:)

第1章矩阵与行列式 >> A(3,:)=A(3,:)+A(2,:) 运行结果： A = 1 2 -1 2 1
第1章矩阵与行列式 >> A(3,:)=A(3,:)+A(2,:) 运行结果： A = >> A(2,:)=-1/3*A(2,:) >> A(1,:)=A(1,:)+A(2,:) 由上可知，方程组有解，其中是自由未知量。

第1章矩阵与行列式解法二：>> A=[2 4 1 1 5;-1 -2 -2 1 -4;1 2 -1 2 1];
第1章矩阵与行列式解法二：>> A=[ ; ; ]; >>A=rref(A) 运行结果： A = 由上可知，结果同解法一。

第1章矩阵与行列式实验四行列式及应用【实验目的】 1. 了解行列式的概念，掌握行列式的性质 2．掌握行列式的计算方法
第1章矩阵与行列式实验四行列式及应用【实验目的】 1. 了解行列式的概念，掌握行列式的性质 2．掌握行列式的计算方法 3．掌握Gramer法则求解线性方程组【实验要求】掌握计算行列式det、解线性方程组solve、生成Vandermonde行列式vander等命令【实验内容】 1．计算下列行列式的值：（1）；（2）；

第1章矩阵与行列式（3）． 2．已知，，验证． 3．用Gramer法则解线性方程组．
第1章矩阵与行列式（3）． 2．已知，，验证． 3．用Gramer法则解线性方程组． 4．c为何值时，齐次线性方程组有非零解？

第1章矩阵与行列式【实验过程】 1．（1）>> A=[ ; ; ; ]; >>det(A) 运行结果： ans = 312 （2）>> A=sym('[a b b b b;b a b b b;b b a b b;b b b a b;b b b b a]'); a^5-10*a^3*b^2+20*a^2*b^3-15*a*b^4+4*b^5 即行列式的值为．（3）>> A=[1,2,3,4]; >> V=vander(A); >> det(V) 12

第1章矩阵与行列式 3．解法一：>> A=[2 1 -5 1;1 4 -7 6;1 -3 0 -6;0 2 -1 2];
第1章矩阵与行列式 3．解法一：>> A=[ ; ; ; ]; >>A1=[ ; ; ; ]; >>A2=[ ; ; ; ]; >>A3=[ ; ; ; ]; >>A4=[ ; ; ; ]; >>a=det(A); >>a1=det(A1);a2=det(A2);a3=det(A3);a4=det(A4); >>X=[a1/a,a2/a,a3/a,a4/a] 运行结果： X = 即得方程组的解为，，，．

第1章矩阵与行列式解法二：>> A=[ ; ; ; ]; >>A1=[A(:,5),A(:,2:4)]; >>A2=[A(:,1),A(:,5),A(:,3:4)]; >>A3=[A(:,1:2),A(:,5),A(:,4)]; >>A4=[A(:,1:3),A(:,5)]; >>a=det(A(:,1:4)); >>a1=det(A1);a2=det(A2);a3=det(A3);a4=det(A4); >>X=[a1/a,a2/a,a3/a,a4/a] 运行结果： X =

第1章矩阵与行列式 4．>> syms c;
第1章矩阵与行列式 4．>> syms c; >>A=[c ;-2 c-1 -2;-2 -2 c-1]; >>a=det(A) ; >> [c]=solve(a,'c') 运行结果： c = [ 5] [ -1] 即当或时，原线性齐次方程组有非零解。

第1章矩阵与行列式实验五向量【实验目的】理解向量、向量的线性组合与线性表示、向量组的线性相关与线性无关的概念
第1章矩阵与行列式实验五向量【实验目的】理解向量、向量的线性组合与线性表示、向量组的线性相关与线性无关的概念掌握向量组线性相关、线性无关的有关性质及判别法理解向量组的极大线性无关组和向量组的秩的概念会求向量组的极大线性无关组和秩 5．掌握矩阵秩的求法【实验要求】掌握简化矩阵为阶梯形式rref、计算行列式det、计算矩阵的秩rank等命令【实验内容】设向量：，，，，问b能否由线性表示？

第1章矩阵与行列式 2．判断下列向量组是否线性相关：（1），，；（2），，． 3．求向量组，，，的一个极大
第1章矩阵与行列式 2．判断下列向量组是否线性相关：（1），，；（2），，． 3．求向量组，，，的一个极大线性无关组，并把其余向量用极大线性无关组中的向量线性表示。 4．求矩阵的秩。

第1章矩阵与行列式 5．求向量在基，，下的坐标．【实验过程】
第1章矩阵与行列式 5．求向量在基，，下的坐标．【实验过程】 1．解法一：>> A=[-1 3 1;0 4 4;1 -2 0;2 5 9]; b=[5;4;-4;1]; B=[A,b]; r=[rank(A),rank(B)] 运行结果： r = 由上可知，故方程组有解。

第1章矩阵与行列式解法二：设，即有 >> A=[-1 3 1 5;0 4 4 4;1 -2 0 -4;2 5 9 1]
第1章矩阵与行列式解法二：设，即有 >> A=[ ; ; ; ] 运行结果： A = >>B=rref(A) B = 由上可知，，故方程组有解，即b可由线性表示，且。

第1章矩阵与行列式 2．（1）>> A1=[1;0;5;7]; >>A2=[-1;1;-2;3];
第1章矩阵与行列式 2．（1）>> A1=[1;0;5;7]; >>A2=[-1;1;-2;3]; >>A3=[2;-1;7;4]; >>A=[A1,A2,A3]; >>r=rank(A) 运行结果： r = 2 此向量组的秩等于2，故此向量组线性相关。（2）>> A1=[1;1;1]; >>A2=[0;2;5]; >>A3=[1;3;6]; >>a=det(A) a = 此向量组组成的矩阵的行列式的值为0，故此向量组线性相关。

第1章矩阵与行列式 3．>> A=[2 3 1 4;1 -1 3 -3;3 2 4 1;-1 0 -2 1];
第1章矩阵与行列式 3．>> A=[ ; ; ; ]; >>B=rref(A) 运行结果： B = 由上可知，，是向量组的一个极大线性无关组，且，． 4．解法一：>> A=[ ; ; ]; >>r=rank(A) r = 2

第1章矩阵与行列式解法二：>> format rat
第1章矩阵与行列式解法二：>> format rat >>A=[ ; ; ]; >>B=rref(A) 运行结果： B = /5 /5 由上可知，矩阵A的秩为2．

第1章矩阵与行列式 5.即求满足方程的解。 >> A1=[1;1;0]; >> A2=[1;0;1];
第1章矩阵与行列式 5.即求满足方程的解。 >> A1=[1;1;0]; >> A2=[1;0;1]; >> A3=[0;1;1]; >> A=[A1,A2,A3]; >> b=[3;-5;9]; >> X=inv(A)*b 输出 X = 8.5000 0.5000

第1章矩阵与行列式设计性实验 Cayler-Hamilton定理【实验目的】 1．理解特征多项式的概念
第1章矩阵与行列式设计性实验 Cayler-Hamilton定理【实验目的】 1．理解特征多项式的概念 2．掌握Cayler-Hamilton定理【实验要求】掌握生成Vandermonde矩阵的vander命令、求矩阵特征多项式系数的poly()命令、求矩阵范数的norm命令及矩阵多项式运算的polyvalm命令【实验内容】 Cayler-Hamilton定理是矩阵理论中的一个比较重要的定理，其内容为：若矩阵A的特征多项式为则有亦即假设矩阵A为Vandermonde矩阵，试验证其满足Cayler-Hamilton定理。

第1章矩阵与行列式【实验方案】 Matlab提供了求取矩阵特征多项式系数的函数poly()，但是poly()函数会产生一定的误差，而该误差在矩阵多项式求解中可能导致了巨大的误差，从而得出错误的结论。在实际应用中还有其他简单的数值方法可以精确地求出矩阵的特征多项式系数。例如，下面给出的Fadeev-Fadeeva递推算法也可以求出矩阵的特征多项式。

第1章矩阵与行列式可以直接由下面的Matlab语句编写一个函数实现： Function c=poly1(A)
第1章矩阵与行列式可以直接由下面的Matlab语句编写一个函数实现： Function c=poly1(A) [nr,nc]=size(A); if nc==nr % 给出若为方阵，则用Fadeev-Fadeeva算法求特征多项式 I=eye(nc); R=I; c=[1 zeros(1,nc)]; for k=1:nc,c(k+1)=-1/k*trace(A*R);r=A*R+c(k+1)*I; end elseif (nr==1 \ nc==1) % 给出为向量时，构造矩阵 A=A(isfinite(A));n=length(A) ; % 出去非数或无界的特征根 c=[1 zeros(1,n)]; for j=1:n c(2:(j+1))=c(2:(j+1))-A(j).*c(1:j); else % 参数有误则给出错误信息 error (’Argument must be a vector or a square matrix.’) end.

第1章矩阵与行列式【实验过程】 >> A = vander([1 2 3 4 5 6 7]); 运行结果： A =
第1章矩阵与行列式【实验过程】 >> A = vander([ ]); 运行结果： A = >> A aa1 = 1.0e+009 *

第1章矩阵与行列式如调用新的poly1()函数，则可以得出如下的精确结果。
第1章矩阵与行列式如调用新的poly1()函数，则可以得出如下的精确结果。 >> aa1=poly1(A);b1=polyvalm(aa1,A);norm(B1) 运行结果： ans = 可见，由此得出的B矩阵就会完全等于0，故该矩阵满足Cayley-Hamilton定理。

第2章线性方程组【线性方程组简介】线性方程组的求解问题促进了线性代数的产生和发展，利用矩阵、行列式和向量这三个基本工具可较好的解决线性方程组的求解问题。利用解向量所构成的基础解系可方便的描述解空间的基本特征及写出通解，从而较好地描述了线性方程组解的结构问题。

第2章线性方程组验证性实验实验一线性方程组【实验目的】理解齐次线性方程组的基础解系、通解及解空间的概念
第2章线性方程组验证性实验实验一线性方程组【实验目的】理解齐次线性方程组的基础解系、通解及解空间的概念掌握齐次线性方程组的基础解系和通解的求法 3．理解非齐次线性方程组解的结构及通解的概念【实验要求】掌握分数数据格式format rat、求基础解系null、简化矩阵为阶梯形式rref、解方程组solve等命令

第2章线性方程组【实验内容】 1.求齐次线性方程组的基础解系及通解。 2.判断方程组是否有解？

第2章线性方程组 3.求方程组的基础解系及通解。 4.求方程组的基础解系及通解。

第2章线性方程组【实验过程】 1．解法一： >> format rat >>A=[ ; ; ] ; >>B=rref(A) 运行结果： B = 1 0 1/2 -1/ /2 3/

第2章线性方程组由上可知，方程组有解，其中，是自由未知量。故得方程组的基础解系为，．通解为，其中为任意常数。

第2章线性方程组解法二： >> format rat >> A=[ ; ; ] ; >> B=null(A,'r') 运行结果： B = -1/2 1/2 -1/2 -3/ >> syms k1 k2 >> X=k1*B(:,1)+k2*B(:,2)

第2章线性方程组运行结果： X = [ -1/2*k1+1/2*k2] [ -1/2*k1-3/2*k2] [ k1] [ k2] 即原方程组的通解为，其中为任意常数。

第2章线性方程组 2．解法一：>> A=[ ; ; ]; >> b=[1;2;4]; >> r=[rank(A),rank([A,b])] 运行结果： r = 2 3 即系数矩阵的秩为2，增广矩阵的秩为3，故原方程组无解。

第2章线性方程组解法二：>> A=[ ; ; ]; >> b=[1;2;4]; >> B=rref([A,b]) 运行结果： B = 从上面所得阶梯形矩阵可以看出系数矩阵的秩为2，增广矩阵的秩为3，故原方程组无解。

第2章线性方程组 3．解法一：>>A=[1 1 1; ; ]; >>b=[66;77;99]; >>r=[rank(A),rank([A,b])] 运行结果： r = 3 3 即系数矩阵的秩等于增广矩阵的秩3，且等于未知量的个数，故原方程组有惟一解。 >> X=inv(A)*b % X=A\b X =

第2章线性方程组解法二：>> syms x1 x2 x3; >>f1=x1+x2+x3-66; >>f2=-10*x1+12*x2+x3-77; >>f3=x1-9*x2+12*x3-99; >> [x1 x2 x3]=solve(f1,f2,f3,'x1','x2','x3') 运行结果： x1 = 21 x2 = 22 x3 = 23

第2章线性方程组 4．解法一：>>A=[ ; ; ]; >>b=[1;2;3]; >>rA=rank(A) 运行结果： rA = 3 >>rAb=rank([A,b]) rAb = 即系数矩阵的秩等于增广矩阵的秩3，故原方程组有解。

第2章线性方程组 >> x0=A\b 运行结果： x0 = 1 0 即原线性方程组的一个特解．

第2章线性方程组 >>B=rref(A) 运行结果： B = 由上可知，原方程组的导出组的解为，即可得其导出组的基础解系为，．故原方程组的通解为，其中为任意常数。

第2章线性方程组解法二： >>A=[1 1 -2 1 3;2 -1 2 2 6;3 2 -4 -3 -9];
第2章线性方程组解法二： >>A=[ ; ; ]; >>b=[1;2;3]; >>X=A\b 运行结果： X = 1

第2章线性方程组 >> B=null(A,'r') 运行结果： B = 故原方程组的通解为，其中为任意常数。

第2章线性方程组设计性实验小行星轨道问题【实验目的】 1. 掌握线性方程组求解 2. 加深对正交变换的理解
第2章线性方程组设计性实验小行星轨道问题【实验目的】 1. 掌握线性方程组求解 2. 加深对正交变换的理解 3. 掌握Matlab软件中的ezplot、zplot命令的区别和适用范围【实验要求】掌握绘制隐函数曲线ezplot命令和彗星状轨迹图comet命令

第2章线性方程组【实验内容】天文学家要确定一颗小行星绕太阳运行的轨道，在轨道平面内建立以太阳为原点的直角坐标系，在两坐标轴上取天文测量单位（一天文单位为地球到太阳的平均距离：9300万里）。在五个不同的时间点对小行星作了观察，测得轨道上五个点的坐标数据如下: 表 2-1 小行星观测数据 x 4.5596 5.0816 5.5546 5.9636 6.2756 y 0.8145 1.3685 1.9895 2.6925 3.5265

第2章线性方程组由开普勒第一定律知，小行星轨道为一椭圆。设方程为试确定椭圆的方程并在轨道的平面内以太阳为原点绘出椭圆曲线。并应用坐标平移变换和正交变换将上例题中的二次曲线方程化为标准方程，绘椭圆轨道图，完成小行星运行的动态模拟。

第2章线性方程组【实验方案】（1）二次曲线方程中有五个待定系数：，，，，。将观察所得的五个点坐标数据，代入二次曲线方程得到关于，，，，的线性方程组求解该方程组得椭圆方程的系数：[ ，，，， ] 。

第2章线性方程组（2）将椭圆的一般方程写成矩阵形式通过变量变换（平移变换和旋转变换）化为椭圆标准方程。首先化去一次项，然后将二次型化为标准型。为了用平移变换消去一次项，令，（，待定），代入方程整理，得

第2章线性方程组其中，。要化简消去一次项，只须选择，使满足二阶线性方程组将，代入椭圆的一般方程，得令求出特征值极其对应的特征向量。可以取与等价的正交单位向量。构造正交矩阵，利用正交变换

第2章线性方程组得椭圆的标准方程：。椭圆长半轴和短半轴分别为，。

第2章线性方程组【实验过程】 (1) MATLAB程序如下：
第2章线性方程组【实验过程】 (1) MATLAB程序如下： x=[4.5596;5.0816;5.5546;5.9636;6.2756]; y=[0.8145;1.3685;1.9895;2.6925;3.5265]; A=[x.^2,2*x.*y,y.^2,2*x,2*y]; b=-[1;1;1;1;1]; a=A\b; syms x y a1 a2 a3 a4 a5 fun=a1*x^2+2*a2*x*y+a3*y^2+2*a4*x+2*a5*y+1; fun=subs(fun,a1,a(1)); fun=subs(fun,a2,a(2)); fun=subs(fun,a3,a(3));

第2章线性方程组 fun=subs(fun,a4,a(4)); fun=subs(fun,a5,a(5)); ezplot(fun,[-1.4,7,-1.5,6.5]) 运行结果： a= 结果表明：二次曲线方程中的各项系数为 = ， =0.1892， = ， =0.4609， =0.4104。

第2章线性方程组图2-2小行星绕太阳运行的轨道

第2章线性方程组 (2) MATLAB程序如下： x=[4.5596;5.0816;5.5546;5.9636;6.2756]; y=[0.8145;1.3685;1.9895;2.6925;3.5265]; A=[x.^2,2*x.*y,y.^2,2*x,2*y]; b=-[1;1;1;1;1]; ak=A\b; C=[ak(1),ak(2);ak(2),ak(3)]; X=-C\[ak(4);ak(5)]; x0=X(1);y0=X(2);X=[X;1]; D=[ak(1),ak(2),ak(4);ak(2),ak(3),ak(5);ak(4),ak(5),1]; F=X'*D*X; [U d]=eig(C);

第2章线性方程组 a=sqrt(-F/d(1,1));b=sqrt(-F/d(2,2)); t=2*pi*(0:5000)/5000; u=a*cos(t);v=b*sin(t); V=U*[u;v]; x1=V(1,:)+x0;y1=V(2,:)+y0; plot(x1,y1,x,y,'*',x0,y0,'rO'),hold on x2=[x1,x1,x1];y2=[y1,y1,y1]; comet(x2,y2) disp([x0,y0]) disp([a,b])

第2章线性方程组图2-3 椭圆轨道图

第2章线性方程组运行结果：。结果表明：椭圆标准方程为：

第3章矩阵的特征值与特征向量【矩阵的特征值与特征向量简介】
第3章矩阵的特征值与特征向量【矩阵的特征值与特征向量简介】矩阵的特征值与特征向量是矩阵的数字特征，利用矩阵的特征值与特征向量可判断矩阵的相似、解决矩阵对角化及实对称矩阵正交化等问题，促进了矩阵理论的进一步发展及应用。

第3章矩阵的特征值与特征向量验证性实验矩阵的特征值与特征向量【实验目的】理解矩阵的特征值和特征向量的概念
第3章矩阵的特征值与特征向量验证性实验矩阵的特征值与特征向量【实验目的】理解矩阵的特征值和特征向量的概念会求矩阵的特征值和特征向量掌握将矩阵化为相似对角矩阵的方法【实验要求】掌握求矩阵的特征多项式poly、求矩阵的特征值和特征向量eig、矩阵的范数norm、值空间正交化orth、单位阵eye等命令

第3章矩阵的特征值与特征向量【实验内容】设，求矩阵A的特征多项式和特征值。求矩阵的特征值和特征向量。
第3章矩阵的特征值与特征向量【实验内容】设，求矩阵A的特征多项式和特征值。求矩阵的特征值和特征向量。求矩阵的特征值和特征向量。

第3章矩阵的特征值与特征向量判断矩阵是否可以对角化，若能，将其对角化。 5．设矩阵，求正交矩阵T，使得为对角矩阵。

第3章矩阵的特征值与特征向量【实验过程】 1．>> >> A=[1 0 0;0 1 8;0 1 3];
第3章矩阵的特征值与特征向量【实验过程】 1．>> >> A=[1 0 0;0 1 8;0 1 3]; >> poly(A) 运行结果： ans = 即矩阵A的特征多项式为． >>lamda=eig(A) lamda = 5 -1 1 即矩阵A的特征值为，，．

第3章矩阵的特征值与特征向量 2．>> A=[1 2 2;2 1 2;2 2 1];
第3章矩阵的特征值与特征向量 2．>> A=[1 2 2;2 1 2;2 2 1]; >> [kesai,lamda]=eig(A) 运行结果： kesai = lamda = 即矩阵A的特征值为，，．对应的特征向量为：，，．

第3章矩阵的特征值与特征向量 3．>> A=sym('[0 0 a;0 b 0;c 0 0]');
第3章矩阵的特征值与特征向量 3．>> A=sym('[0 0 a;0 b 0;c 0 0]'); >> [kesai,lamda]=eig(A) 运行结果： kesai = [ 1/c*(c*a)^(1/2), -1/c*(c*a)^(1/2), ] [ , , ] [ , , ] lamda = [ (c*a)^(1/2), , ] [ , -(c*a)^(1/2), ] [ , , b] 即矩阵A的特征值为，，．对应的特征向量为：，，．

第3章矩阵的特征值与特征向量 4．>> A=[0 1 1;1 0 1;1 1 0];
第3章矩阵的特征值与特征向量 4．>> A=[0 1 1;1 0 1;1 1 0]; >> [kesai,lamda]=eig(A) 运行结果： kesai = lamda = >> inv(kesai)*A*kesai-lamda ans = 1.0e-015 *

第3章矩阵的特征值与特征向量因1.0e-015<<1，故A能相似于对角矩阵。 >> kesai'*A*kesai
第3章矩阵的特征值与特征向量因1.0e-015<<1，故A能相似于对角矩阵。 >> kesai'*A*kesai 运行结果： ans = 即存在正交矩阵使得

第3章矩阵的特征值与特征向量 5．解法一：>> A=[ ; ; ; ] ; >> [kesai,lamda]=eig(A) 运行结果： kesai = lamda = 即所求正交矩阵为．

第3章矩阵的特征值与特征向量 >> kesai'*A*kesai 运行结果： ans = 4 * * * * 4 * *
第3章矩阵的特征值与特征向量 >> kesai'*A*kesai 运行结果： ans = * * * * * * * * * * 即经验证有． >> norm(kesai'*kesai-eye(4)) 9.7171e-016 由上可知，所求正交矩阵精度很高。

第3章矩阵的特征值与特征向量解法二：>> A=[ ; ; ; ] ; >>T=orth(A) 运行结果： T = >> norm(T'*T-eye(4)) ans = 7.6679e-016 由上可知，所求正交矩阵精度很高。

第3章矩阵的特征值与特征向量设计性实验【实验目的】 1.掌握矩阵的相似变换
第3章矩阵的特征值与特征向量设计性实验实验一矩阵相似变换在控制理论中的应用【实验目的】 1.掌握矩阵的相似变换 2.利用矩阵相似变换方法，将控制理论中一般的状态方程变换成某种特殊的形式，以便于更好地进行系统的性质分析 3.掌握控制系统的可控标准型、可观察标准型和Jordan标准型【实验要求】掌握Matlab软件中有关相似变换的命令【实验内容】给出系统的相似变换的概念，介绍基于矩阵相似变换的各种标准及变换方法，并用MATLAB编程实现。试求出下面系统的可控标准型: ，并求该状态方程模型的可观测标准型以及Jordan标准型。

第3章矩阵的特征值与特征向量【实验方案】 1.线性系统的相似变换
第3章矩阵的特征值与特征向量【实验方案】 1.线性系统的相似变换假设存在一个非奇异矩阵，且定义了一个新的状态变量使得，这样关于新状态变量的状态方程模型可以写成，且式中，，，。在矩阵下的状态变换称为相似变换，称为相似变换矩阵。

第3章矩阵的特征值与特征向量 2.单变量控制系统的可控、可观测标准型转换对单变量系统（1）来说，若系统的特征多项式可以写成 +
第3章矩阵的特征值与特征向量 2.单变量控制系统的可控、可观测标准型转换对单变量系统（1）来说，若系统的特征多项式可以写成 + 则可以够造出变换矩阵这样就可以将原来系统变换成可控制标准型。可以用容易地写出变换矩阵 ; %求特征多项式系数，建议用ploy1()取代ploy()

第3章矩阵的特征值与特征向量 3.控制系统的Jordan标准型转换
第3章矩阵的特征值与特征向量 3.控制系统的Jordan标准型转换系统的Jordan标准型可以由函数直接求出。值得指出的是，若系统的矩阵含有复数特征值，则用函数不能得出正确结果，应该结合前面Jordan变换的方法手工构造变换矩阵，得出合适的变换系统。

第3章矩阵的特征值与特征向量实验二矩阵的三角分解【实验目的】 1.理解矩阵的三角分解（又称为LU分解） 2.掌握函数的两种调用方法【实验要求】掌握Matlab软件中有关矩阵LU分解的命令【实验内容】分别用两种方法调用MATLAB中的函数，实现矩阵LU分解问题。

第3章矩阵的特征值与特征向量【实验方案】矩阵的三角分解又称为LU分解，它的目的是将一个矩阵分解成一个下三角矩阵L和一个上三角矩阵U的乘积，亦即A=LU,其中L和U矩阵可以分别写成由这两个矩阵可以简单的写出一个矩阵 ,其中

第3章矩阵的特征值与特征向量这样产生的矩阵与原来的A矩阵的关系可以写成

第3章矩阵的特征值与特征向量因此，可以立即得出求取和的递推计算公式该公式的递推初值为注意，在上述的算法中并未对主元素进行任何选取，因此该算法并不一定数值稳定，因为在运算过程中0或很小的数值可能被用作除数。在MATLAB中也给出了基于主元素的矩阵LU分解函数，该函数的调用格式为

第3章矩阵的特征值与特征向量分解，为置换矩阵，其中，分别为变换后的下三角和上三角矩阵。在MATLAB的函数中考虑了主元素选取的问题，所以该函数一般会给出可靠的结果。又该函数得出的下三角矩阵L并不一定是一个真正的下三角矩阵，因为选取它可能进行了一些元素行的交换，这样主对角线的元素可能不是1，而在矩阵L内存在一个唯一的置换，其各个元素的值均是1.如果想获得有关换行信息，则可以由后一种格式调用函数，这时P为单位阵变换出的置换矩阵，A矩阵可以分解成。

第3章矩阵的特征值与特征向量【实验过程】（1）求出三角分解矩阵。

第3章矩阵的特征值与特征向量可见，这样得出的矩阵并非下三角矩阵，这是因为再分解过程中采用了主元素交换的方法。现在考虑函数的另一中调用方法。

第3章矩阵的特征值与特征向量注意，这里得出的P矩阵不是一个单位矩阵，而是单位矩阵的置换矩阵。结合得出的矩阵可以看出，P矩阵的，表明需要将矩阵的第4行换到第2行，表明需要将的第2行换至第3行，将原来第3行换至第4行，这样就可以得出一个真正的下三角矩阵L了。将L,P,U代入并检验，可以精确地还原A矩阵。

第3章矩阵的特征值与特征向量

第3章矩阵的特征值与特征向量实验三线性方程组的数值解法【实验目的】 1．掌握求解上三角形方程组的回代法，求解线性方程组的Gauss消去法、LU分解、Jacobi迭代法、Gauss-Seidel迭代法 2．比较各种算法的收敛条件及运算效率 3．使用Matlab编写线性方程组数值解法的M文件及程序【实验要求】 1．求解上三角形方程组的回代法、求解线性方程组的Gauss消去法、LU分解、Jacobi迭代法、Gauss-Seidel迭代法的原理及步骤 2．能Matlab编写线性方程组数值解法的M文件及程序

第3章矩阵的特征值与特征向量【实验内容】用回代法解上三角形方程组，用列主元Gauss消去法、矩阵直接三角分解法、Jacobi迭代法和Gauss-Seidel迭代法解线性方程组。 1.求解上三角线性方程组 2.用列主元Gauss消去法求解线性方程组

第3章矩阵的特征值与特征向量 3.用矩阵的LU 分解求解线性方程组 4.用Jacobi迭代法求解线性方程组 5.用Gauss-Seidel迭代法求解线性方程组

第3章矩阵的特征值与特征向量【实验方案】对于上三角形方程组，其矩阵形式表示为：形如：其中，
第3章矩阵的特征值与特征向量【实验方案】对于上三角形方程组，其矩阵形式表示为：形如：其中， U称为上三角矩阵。若，即，亦即矩阵 U是非奇异的，则原方程组有唯一解，且可从最后一个方程解出，即：

第3章矩阵的特征值与特征向量代入倒数第二个方程，得到：一般地，设已求得 ,则由方程的第 i个方程，可得：
第3章矩阵的特征值与特征向量代入倒数第二个方程，得到：一般地，设已求得 ,则由方程的第 i个方程，可得：上述求解过种，称为回代过程。回代过程所用乘除法运算次数为，加减法运算次数为

第3章矩阵的特征值与特征向量 2. Gauss消去法的操作过程如下：将原线性方程组写成增广矩阵的形式，其中， , , 第一步消元：设，将中的第一行乘以，加到第 i行上去，，可得到同解方程组的增广矩阵：

第3章矩阵的特征值与特征向量其中，通常，称为消元因子，表为主元。上述做法，直至第 n-1步完成，得到同解的方程组：

第3章矩阵的特征值与特征向量即为上三角矩阵，可表示为式(11.2.2)，这时用回代过程便可求得解向量 x。而列主元Gauss消去法的基本思想是：设已作 k-1步消元，在进行第 k步消元之前，选出第 k列中位于对角线及其以下元素绝对值中的最大者，即确定t ，使得：将的第t行和第k行互相交换，则元素为新的主元素，其余元素也均以交换后的位置表示。然后再按Gauss消去法进行第k 步消元。这种方法称为列主元Gauss消去法。此时，消元因子满足：

第3章矩阵的特征值与特征向量一般能保证舍入误差不扩散，这个方法基本是稳定的。列主元Gauss消去法的运算量除选主元及行交换外和Gauss消去法是相同的。

第3章矩阵的特征值与特征向量 3. Gauss消去过程，实际就是对方程组的增广矩阵连续左乘以得到上三角方程组，即，其中：
第3章矩阵的特征值与特征向量 3. Gauss消去过程，实际就是对方程组的增广矩阵连续左乘以得到上三角方程组，即，其中：令，用左乘(11.3.1)式，得，

第3章矩阵的特征值与特征向量于是，有：称为单位下三角矩阵，为上三角矩阵。这种矩阵的分解为简单矩阵的乘积，称为矩阵的分解。

第3章矩阵的特征值与特征向量 4. 对于线性方程组设。从第个方程组解出，得到如下同解方程组：

第3章矩阵的特征值与特征向量建立相应迭代格式：称上式为Jacobi迭代格式。

第3章矩阵的特征值与特征向量 5. Jacobi迭代格式在逐个求的分量时，当计算到时，分量都已求出，但却被束之高阁，而仍用旧分量进行计算。直观上看，最新算出的分量可能比旧的分量要准确些。因此，设想一旦当新分量已求出，马上就用它来替代，也就是在Jacobi迭代法中求时用代替，这就是Gauss-Seidel迭代法。

第3章矩阵的特征值与特征向量【实验过程】建立回代法的M函数文件如下： function X=backsub(A,b)
第3章矩阵的特征值与特征向量【实验过程】建立回代法的M函数文件如下： function X=backsub(A,b) % A是一个n阶上三角非奇异阵。 % b是一个n维向量。 % X是线性方程组AX=b的解。 n=length(b); X=zeros(n,1); X(n)=b(n)/A(n,n); for k=n-1:-1:1 X(k)=(b(k)-A(k,k+1:n)*X(k+1:n))/A(k,k); end 然后，用上述程序求解上三角线性方程组：

第3章矩阵的特征值与特征向量建立一个主程序prog1121.m clc clear
第3章矩阵的特征值与特征向量建立一个主程序prog1121.m clc clear A=[3,-2,1,-4;0,4,-1,2;0,0,2,3;0,0,0,5]; b=[8;-3;11;15]; backsub(A,b) 然后在MATLAB命令窗口运行上述主程序，即： >> prog1121

第3章矩阵的特征值与特征向量运行结果： ans = 5 -2 1 3 建立列主元 Gauss消去法的函数文件如下：
第3章矩阵的特征值与特征向量运行结果： ans = 5 -2 1 3 建立列主元 Gauss消去法的函数文件如下： function X=uptrbk(A,b) % A是一个n阶矩阵。 % b是一个n维向量。 % X是线性方程组AX=b的解。

第3章矩阵的特征值与特征向量 [N N]=size(A); X=zeros(1,N+1); Aug=[A b]; for p=1:N-1
第3章矩阵的特征值与特征向量 [N N]=size(A); X=zeros(1,N+1); Aug=[A b]; for p=1:N-1 [Y,j]=max(abs(Aug(p:N,p))); C=Aug(p,:); Aug(p,:)=Aug(j+p-1,:); Aug(j+p-1,:)=C; if Aug(p,p)==0 'A是奇异阵，方程无惟一解' break end

第3章矩阵的特征值与特征向量 for k=p+1:N m=Aug(k,p)/Aug(p,p);
第3章矩阵的特征值与特征向量 for k=p+1:N m=Aug(k,p)/Aug(p,p); Aug(k,p:N+1)=Aug(k,p:N+1)-m*Aug(p,p:N+1); end % 这里用前面程序中定义的函数backsub来进行回代。 X=backsub(Aug(1:N,1:N)，Aug(1:N,N+1)；

第3章矩阵的特征值与特征向量再用上述程序求解线性方程组：建立一个主程序prog1122.m clc clear
第3章矩阵的特征值与特征向量再用上述程序求解线性方程组：建立一个主程序prog1122.m clc clear A=[2,4,-6;1,5,3;1,3,2]; b=[-4;10;5]; uptrbk(A,b) 然后在MATLAB命令窗口运行上述主程序，即：

第3章矩阵的特征值与特征向量 >> prog1122 运行结果： ans = -3 2 1
第3章矩阵的特征值与特征向量 >> prog1122 运行结果： ans = -3 2 1 建立矩阵分解解方程组的函数文件如下： function X=lufact(A,b) % A为n阶矩阵。 % b是n维向量。 % X是所求的AX=b的解。

第3章矩阵的特征值与特征向量 [N,N]=size(A); X=zeros(N,1); Y=zeros(N,1);
第3章矩阵的特征值与特征向量 [N,N]=size(A); X=zeros(N,1); Y=zeros(N,1); C=zeros(1,N); R=1:N; for p=1:N-1 [max1,j]=max(abs(A(p:N,p))); C=A(p,:); A(p,:)=A(j+p-1,:); A(j+p-1,:)=C; d=R(p); R(p)=R(j+p-1); R(j+p-1)=d;

第3章矩阵的特征值与特征向量 if A(p,p)==0 'A是奇异阵，方程组无惟一解' break end for k=p+1:N
第3章矩阵的特征值与特征向量 if A(p,p)==0 'A是奇异阵，方程组无惟一解' break end for k=p+1:N mult=A(k,p)/A(p,p); A(k,p)=mult; A(k,p+1:N)=A(k,p+1:N)-mult*A(p,p+1:N); Y(1)=b(R(1)); for k=2:N Y(k)=b(R(k))-A(k,1:k-1)*Y(1:k-1);

第3章矩阵的特征值与特征向量 X(N)=Y(N)/A(N,N); for k=N-1:-1:1
第3章矩阵的特征值与特征向量 X(N)=Y(N)/A(N,N); for k=N-1:-1:1 X(k)=(Y(k)-A(k,k+1:N)*X(k+1:N))/A(k,k); end 用上述程序求解线性方程组：建立一个主程序prog1131.m clc clear A=[1,3,5,7;2,-1,3,5;0,0,2,5;-2,-6,-3,1]; b=[1,2,3,4]; Lufact(A,b)

第3章矩阵的特征值与特征向量 >> prog1131 运行结果： ans = 1.3429 0.6857.0000
第3章矩阵的特征值与特征向量然后在MATLAB命令窗口运行上述主程序，即： >> prog1131 运行结果： ans = 1.3429 1.8000 -3

第3章矩阵的特征值与特征向量还有另一种简便的解法是直接使用MATLAB中对矩阵进行LU分解的命令，然后再用回代法求解原方程组，过程如下：
第3章矩阵的特征值与特征向量还有另一种简便的解法是直接使用MATLAB中对矩阵进行LU分解的命令，然后再用回代法求解原方程组，过程如下： clc clear A=[1,3,5,7;2,-1,3,5;0,0,2,5;-2,-6,-3,1]; b=[1,2,3,4]; %对系数矩阵A进行LU分解 [L,U]=lu(A); %用前面程序中定义的函数backsub来进行回代 backsub(U,b)

第3章矩阵的特征值与特征向量建立Jacobi迭代格式的函数文件如下：
第3章矩阵的特征值与特征向量建立Jacobi迭代格式的函数文件如下： function X=jacobi(A,b,P,delta,max1) % A是n维非奇异阵。 % b是n维向量。 % P是初值。 % delta是误差界。 % max1是给定的迭代最高次数。 %X为所求的方程组AX=b的近似解。

第3章矩阵的特征值与特征向量 N=length(b); for k=1:max1 for j=1:N
第3章矩阵的特征值与特征向量 N=length(b); for k=1:max1 for j=1:N X(j)=(b(j)-A(j,[1:j-1,j+1:N])*P([1:j-1,j+1:N]))/A(j,j); end err=abs(norm(X'-P)); P=X'； if(err<delta) break X=X';k,err;

第3章矩阵的特征值与特征向量用上述程序求解线性方程组：建立一个主程序prog1141.m clc clear
第3章矩阵的特征值与特征向量用上述程序求解线性方程组：建立一个主程序prog1141.m clc clear A=[4,1,-1;1,-5,-1;2,-1,-6]; b=[13;-8;-2]; P=[0;0;0]; jacobi(A,b,P,10^(-4),20)

第3章矩阵的特征值与特征向量然后在MATLAB命令窗口运行上述主程序，即： >> prog1141 运行结果： k = 9
第3章矩阵的特征值与特征向量然后在MATLAB命令窗口运行上述主程序，即： >> prog1141 运行结果： k = 9 ans = 3.0000 2.0000 1.0000

第3章矩阵的特征值与特征向量建立Gauss-Seidel迭代法的M函数文件如下：
第3章矩阵的特征值与特征向量建立Gauss-Seidel迭代法的M函数文件如下： function X=gseid(A,b,P,delta,max1) % A是n维非奇异阵。 % b是n维向量。 % P是初值。 % delta是误差界。 % max1是给定的迭代最高次数。 % X为所求的方程组AX=b的近似解。

第3章矩阵的特征值与特征向量 N=length(b); for k=1:max1 for j=1:N if j==1
第3章矩阵的特征值与特征向量 N=length(b); for k=1:max1 for j=1:N if j==1 X(1)=(b(1)-A(1,2:N)*P(2:N))/A(1,1); elseif j==N X(N)=(b(N)-A(N,1:N-1)*(X(1:N-1))')/A(N,N); else X(j)=(b(j)-A(j,1:j-1)*X(1:j-1)-A(j,j+1:N)*P(j+1:N))/A(j,j); end err=abs(norm(X'-P)); relerr=err/(norm(X)+eps); P=X';

第3章矩阵的特征值与特征向量 if(err<delta)|(relerr<delta) break end X=X‘;
第3章矩阵的特征值与特征向量 if(err<delta)|(relerr<delta) break end X=X‘; err,k 然后，用上述程序求解线性方程组：建立一个主程序prog1142.m

第3章矩阵的特征值与特征向量 clear A=[4,-1,1; 4,-8,1; -2,1,5] b=[7;-21;15]
第3章矩阵的特征值与特征向量 clear A=[4,-1,1; 4,-8,1; -2,1,5] b=[7;-21;15] P=[1;2;2] gseid(A,b,P,10^(-4),14) 然后在MATLAB命令窗口运行上述主程序，即： >> prog1142 运行结果： err = e-005 k = 6 ans = 2.0000 4.0000 3.0000 经过6次迭代得到满足精度要求的解。

第3章矩阵的特征值与特征向量实验四房屋装修的工资问题【实验目的】 1．理解矩阵特征值概念
第3章矩阵的特征值与特征向量实验四房屋装修的工资问题【实验目的】 1．理解矩阵特征值概念 2．能根据实际问题，建立模型然后使用Matlab相关命令求解【实验要求】掌握求解特征值的eig命令、生成对角矩阵的diag命令等

第3章矩阵的特征值与特征向量【实验内容】有三个技术个人分别是木工、电工和管道工，他们准备合作装修自己的新房子。在装修之前约定：每人总共工作20天（包括在自己家）；每人每日的工资平均为100元；每人的日工资应使得每人的总收入和总支出等。需要计算每人的日工资分别是多少，以确定他们的工作日交换是否平衡，如果不平衡，将由谁买单。一个初步的工作日分配方案如下表3-1 工作日分配方案工作日工种木工电工管道工木工家 4 2 12 电工家 8 10 管道工家 6

第3章矩阵的特征值与特征向量【实验方案】设木工、电工和管道工的日工资分别为：，，。由总收入和总支出相等的约定，建立线性议程组整理，得

第3章矩阵的特征值与特征向量显然问题与矩阵特征值问题有联系，由于矩阵
第3章矩阵的特征值与特征向量显然问题与矩阵特征值问题有联系，由于矩阵是正矩阵且每列元素之和均为20，所以20是该矩阵的牲值，于是就是属于特征值的特征向量。按约定总工作量决定总工资应该为6000元，则应该有

第3章矩阵的特征值与特征向量【实验过程】 MATLAB程序如下 A=[4,2,12;8,10,2;8,8,6];
第3章矩阵的特征值与特征向量【实验过程】 MATLAB程序如下 A=[4,2,12;8,10,2;8,8,6]; [P,D]=eig(A); disp(diag(D)) II=input('input Index about eigvalu=20:='); if II==0,error('problem have no solution'),end alpha=P(:,II); R=alpha./sum(alpha); format bank daily=300*R' pay=A*diag(daily)

第3章矩阵的特征值与特征向量运行结果：在MATLAB命令窗口中运行程序，屏幕将显示出A的三个特征值 20.00 -2.00 2.00
第3章矩阵的特征值与特征向量运行结果：在MATLAB命令窗口中运行程序，屏幕将显示出A的三个特征值 20.00 -2.00 2.00 由于第一个特征值恰好为20，在提示符“input Index about eigvalu=20:=”后输入索引值1。程序继续运行，得出最后计算结果为 daily = pay = 每人的日工资由变量daily的数据给出。

第3章矩阵的特征值与特征向量结果表明：表3-2 日工资列表最后的二维数组给出了二维数组，表明付款明细账，行表示支付，列表示收取。显然第一行相加等于第一列相加，第二行相加等于第二列相加，第三行相加等于第三列相加。表3-3 工资支付收取方案工种木工电工管道工日工资 93.94 96.97 109.09 支付收取木工电工管道工 375.76 193.94 751.52 969.70 218.18 775.76 654.55

第3章矩阵的特征值与特征向量实验五层次分析法【实验目的】 1、掌握正互反矩阵、一致性矩阵等一些基本知识
第3章矩阵的特征值与特征向量实验五层次分析法【实验目的】 1、掌握正互反矩阵、一致性矩阵等一些基本知识 2、理解层次分析法、理想解逼近排序法的基本原理和步骤【实验要求】 1、能够对一些定性问题或半定性问题能用层次分析法进行建模求解,并用MATLAB编程实现 2、掌握用层次分析法进行建模求解过程Matlab中一些相关命令的用法

第3章矩阵的特征值与特征向量【实验内容】中国石化股份有限公司西南分公司日前开采的气田主要有孝泉气田、新场气田、合兴场气田等，这些气田最早于1988年投入开采，现已形成较为完善的集输管网系统，日输气能力已达到近400万方/天。川西气田天然气气质好，甲烷含量高达90％以上，不含H2S，因此从气井产出的天然气一般只经过采集气站分离、计量、调压后直接外输进入集输管网供给用户。但随着气田开发的进行，气井压力降低，产水逐渐增加，越来越多的气井采用泡沫排水采气工艺进行排水采气，由于站内流程分离器分离效果不佳，加之气井采用泡沫排水采气，部分水分形成泡沫混合在天然气中，导致分离器的分离效果降低，导致天然气含水量增大，对整个集输管网系统运行带来了一系列不利影响。因此，开展川西气田集输管网系统脱水技术应用研究，对提高用户供气质量，降低管网运行压损，消除管内积液对管道的腐蚀具有极为重要的作用。根据不同脱水装置设置点、选

第3章矩阵的特征值与特征向量址方案以及不同的脱水方法（三甘醇脱水法和硅胶脱水法），得到了6种川西气田集输管网整体脱水方案，见表3-4。
第3章矩阵的特征值与特征向量址方案以及不同的脱水方法（三甘醇脱水法和硅胶脱水法），得到了6种川西气田集输管网整体脱水方案，见表3-4。表3-4 各个方案的特征指标表指标方案总投资 (万元) 运营成本 (万元/年) 处理量（万m3/日）压损（MPa）露点（℃）保护管线长度（公里）方案1 261.90 435 0.2 -5 方案2 256.41 0.4 -10 方案3 222.96 326.4 398 方案4 218.64 方案5 222.43 400 379 方案6 183.98

第3章矩阵的特征值与特征向量由于6个方案都有各自的优点和缺点，如果单纯从技术或经济某一个方面来决定选择一个方案作为最终方案显得不太合理，因此，需要对6个方案进行技术经济综合评价，选择技术经济综合评价相对较高的方案作为最终方案。多方案的技术经济综合评价的主要思路：首先选取综合评价特征指标，由于各特征指标重要性并不完全相同，因而通过层次分析法把各指标重要性从定性转化成定量，作为权重；再通过逼近理想解排序法（TOPSIS法）将各个方案排出优劣次序，确定推荐方案。

第3章矩阵的特征值与特征向量【实验方案】层次分析法的基本原理和步骤：层次分析法的基本思路是：将研究对象分解为不同的组成因素,按各因素之间的隶属关系,把它们排成从高到低的若干层次,建立递阶层次结构。对同层的各元素进行两两比较,就每一层次的相对重要性予以定量表示,并利用数学方法确定出每一层次各项因素的权值。层次分析法的流程图如下图所示，步骤的详细描述请参见文献《数学模型（第三版）》姜启源，谢金星。

第3章矩阵的特征值与特征向量 N Y 建立层次结构模型构造判断矩阵计算特征值和特征向量一致性检验通过吗？确定指标权重
第3章矩阵的特征值与特征向量建立层次结构模型构造判断矩阵计算特征值和特征向量一致性检验通过吗？确定指标权重 N Y 三标度层次分析法流程图

第3章矩阵的特征值与特征向量理想解逼近排序法（TOPSIS法）的基本原理和步骤：理想解逼近排序法的基本原理为：在维空间中，将方案集中的各个备选方案与理想方案和负理想方案的距离进行比较，既靠近理想方案又远离负理想方案的方案是方案集中的最佳方案；并可以据此排定方案集中各备选方案的优先序。 TOPSIS法的流程图如下图所示，步骤的详细描述请参阅相关文献。

第3章矩阵的特征值与特征向量 1.建立层次分析法模型 (1) 建立层次结构模型：根据专家意见，选定方案总投资、年运营成本、方案总脱水处理量、脱水前后压损、脱水后干气露点值、保护管线总长为综合评价指标，采用三标度层次分析法构成层次分析模型，对方案进行技术经济评价。综合评价层次结构图3-7：

第3章矩阵的特征值与特征向量 (2)构造判断矩阵：从层次结构模型的第2层开始，对于从属于（或影响及）上一层每个因素的同一层诸因素，用成对比较法构造判断矩阵。在进行成对比较时，可采用Saaty等人提出的1-9比较尺度。但这里采用了一种较为简化的构造方法：确定各指标综合按重要性大小排列顺序如下：总投资> 处理量>运营成本> 压损>露点>保护管线长度，即按照层次结构图特征指标顺序，设为第 i个特征指标（ i=1,…,6），构造出判断矩阵为（其中特征指标比大则为2，相等则为1，小则为0.5），即：

第3章矩阵的特征值与特征向量 (3)求特征值和特征向量：直接通过计算机编程求解矩阵的特征值和特征向量，找出最大的特征值和其对应的特征向量，最大特征值为，其对应的特征向量为： (0.6283，0.3958，0.4986，0.3141，0.2493，0.1979)T 计算矩阵特征值和特征向量的matlab命令： >> A=[1,2,2,2,2,2; 0.5,1,0.5,2,2,2; 0.5,2,1,2,2,2; 0.5,0.5,0.5,1,2,2; 0.5,0.5,0.5,0.5,1,2; 0.5,0.5,0.5,0.5,0.5,1]; >> [V,D]=eig(A); >>B=V(:,1)

第3章矩阵的特征值与特征向量 (4)一致性检验：一致性指标的计算公式为：，其中为最大特征根，为矩阵的阶。计算得出，故。根据由Saaty等人给出的层次分析一致性表查得，维数为6的随机一致性指标为1.24，所以，检验通过，认为判断矩阵一致性可以接受。其权重为对应向量，将权重化为标准权重，即为：0.2751，0.1733，0.2183，0.1375，0.1092，可得到标准权重如表3-8:

第3章矩阵的特征值与特征向量表3-8 特征指标的标准权重表计算标准权重的matlab命令： >> B=B/sum(B)

第3章矩阵的特征值与特征向量 2 理想解逼近排序法（TOPSIS法）（1）计算规范决策矩阵：首先对原始数据（表3-1）标准化处理，得到无量纲数据；然后用向量规范化的方法求得规范决策矩阵。标准化处理采用：经过以上公式标准化处理后得到矩阵

第3章矩阵的特征值与特征向量在matlab中编写标准化处理的m文件：

第3章矩阵的特征值与特征向量 clc clear y=[ ; ; ; ; ; ]; c=zeros(6,6); for j=1:6 for i=1:6 c(i,j)=y(i,j)/sqrt(sum(y(:,j).^2)); end c

第3章矩阵的特征值与特征向量（2）构成加权规范阵。设由层次分析法得到的权重为则，其中，即得到为：

第3章矩阵的特征值与特征向量计算的matlab程序命令： for i=1:6 for j=1:6 x(i,j)=w(1,j)*c(i,j); end

第3章矩阵的特征值与特征向量 (3)确定理想方案和负理想方案第个属性值为，则 j为效益型属性（值越大越好）理想解 j为成本型属性（值越小越好）负理想解理想方案的特征值分别为： =（0.0878,0.0568,0.0995,0.0355, ,0.0444)T 负理想方案的特征值分别为： = (0.1279, , , , , 0) T

第3章矩阵的特征值与特征向量（4）计算各个方案到理想方案和到负理想方案的距离。备选方案到理想方案的距离为得到6个方案到理想方案的距离为：0.5089，0.5397，0.5134，0.5664，0.5276，个方案到负理想方案的距离为：0.4550，0.5211，0.5382，0.5460，0.5410，

第3章矩阵的特征值与特征向量得到6个方案到负理想方案的距离为：计算距离的matlab命令： x1=[0.0878,0.0568,0.0995,0.0355, ,0.0444]; for i=1:6 d1(i)=sqrt(sum(x(i,:)-x1(i).^2)); end x2=[0.1279, , , , , 0]; d2(i)=sqrt(sum(x(i,:)-x2(i).^2));

第3章矩阵的特征值与特征向量（5）计算各方案的排序指标值（即综合评价指数）：各方案的综合指标、、、、、分别为：，0.9654，1.0482，0.9640，1.0254，（6）按由大到小排序如下： > > > > > 根据其原理，综合指标越大，表示到差理想方案与到理想方案的距离之比越大，故方案越好。因此，方案好坏的次序为： > > > > >

第3章矩阵的特征值与特征向量其综合评价图形如下：图3-9综合评价方案直观图

第3章矩阵的特征值与特征向量在图3-9中，横坐标表示备选方案到理想方案的距离，纵坐标表示到负理想方案的距离，直线与横坐标的夹角的正切为综合指标，其夹角越大表示方案越好。因此，方案3最佳。

第3章矩阵的特征值与特征向量实验六卷烟叶组配方设计【实验目的】掌握线性方程组的各种解法。能根据实际问题，使用Matlab建立相应的线性方程组并求解。【实验要求】 1．掌握几种线性方程组（定解方程组、不定方程组、超定方程组、奇异方程组、符号方程组）的解法。 2．能用Matlab求解不同类型线性方程组的方法。

第3章矩阵的特征值与特征向量【实验内容】如何提高卷烟抽吸时的感官质量，以及如何降低烟气中的有害成分始终是卷烟制造工业的重中之中。卷烟的叶组配方，即卷烟中的混合烟丝，是由多种单料烟叶按照某种特定的百分比例组合而成的，其化学成分含量（包括总糖、总碱、氯、氮、磷、氧化钾的含量）与其感官质量指标（包括光泽、香气、谐调、杂气、刺激性、余味）和烟气化学成分含量（包括焦油量、CO量、烟气烟碱量）之间存在着一定的映射关系，也就是说特定化学成分的叶组配方对应着其特定的感官质量和烟气化学成分，叶组配方化学成分的含量从另一个角度反映了其感官质量和烟气化学成分含量。因此，在对叶组配方进行设计时，通常要求在确定叶组配方化学成分含量的前提下来确定进入叶组的各

第3章矩阵的特征值与特征向量种单料烟叶的百分比例，这样既保证了卷烟叶组的感官质量，又确保了其烟气化学成分含量不会太高。本实验要求设计出根据叶组配方化学成分含量要求确定各种单料烟叶百分比例的数学模型，并对模型求解，给出问题的结果。配方设计师根据将要生产的卷烟的抽吸风格的需要选择了12种单料烟叶进入叶组，其化学成分含量与叶组所要求的化学成分含量如表3-10所示，要求根据叶组所要求的化学成分含量确定出各种单料烟叶在叶组中所含的百分比例。

第3章矩阵的特征值与特征向量【实验方案】该问题的目的是要在确定叶组化学成分的前提下求出每种单料烟叶在叶组中所占的百分比例，而叶组的某种化学成分是由各种单料烟叶相对应的化学成分按照其百分比例组合而成的，并且各种单料烟叶的百分比例之和应该为100%，因此，我们可据此列出线性方程组，求解该线性方程组即可求得每种单料烟叶的百分比例。设编号为i的单料烟叶在该叶组中所占的百分比例为，即编号依次为1，2，…，12的单料烟叶在叶组中所占的百分比例分别为，，…，，根据前面表中给出的数据可列出下面的线性方程组

第3章矩阵的特征值与特征向量求解该线性方程组即可求得12种单料烟叶的百分比例。
第3章矩阵的特征值与特征向量求解该线性方程组即可求得12种单料烟叶的百分比例。观察前面所列出的方程组，未知数的个数大于方程组的个数，该线性方程组是不定方程组，有多个解，可利用线性代数中求解不定线性方程组的方法，求出该方程组的特解与通解。

第3章矩阵的特征值与特征向量【实验过程】 >> clear all >> clc >> %输入方程组的系数矩阵 >>A=[31.75,32.49,26.64,25.5,28.31,24.13,22.67,18.72,19.48,24.76,19.36,22.44; 2.96,2.07,1.77, 2.35,1.76,2.91,2.59,3,2.07,2.47,3.41,2.48; 0.24,0.21,0.22,0.22,0.3 ,0.13,0.13,0.07,0.36,0.27,0.26,0.45; 2.05,1.84,1.86,2.14,1.99,2.36,1.97,2.22,1.79,2.06,2.21,2.39; 0.2,0.22,0.26,0.25,0.27,0.2,0.21,0.2,0.26,0.21,0.2,0.19; 2.01,1.95,2.29,2.11,2.1,2.36,2.37,2.25,1.89,2.26,2.22,2.43; 1,1,1,1,1,1,1,1,1,1,1,1];

第3章矩阵的特征值与特征向量 >> B=[25.26;2.43;0.23;2.05;0.23;2.18;1]; >> X0=A\B %求方程组的一个特解运行结果： X0 =

第3章矩阵的特征值与特征向量 >> X=null(A) %求原方程组对应的齐次线性方程组的基础解系运行结果： X =

第3章矩阵的特征值与特征向量所以，原方程组的通解为

第3章矩阵的特征值与特征向量如果要求出原方程组的一个特定解，即求出叶组配方中每种单料烟叶的具体百分比例，就将 , , , , 取某特定的值，代入上式进行计算就可以了，但一定要确保计算出的特定解的所有元素都要大于0且小于1，因为每种单料烟叶在叶组中所占的百分比例必定是大于0且小于1的。如取 =-0.001， =0.15， =0.001， =-0.01， =-0.2，计算出的结果为x1=0.0143，x2=0.1241，x3=0.0952，x4=0.0670，x5=0.2308，x6=0.0811，x7=0.1010，x8=0.0311，x9=0.0165，x10=0.0370，x11=0.1952，x12=0.0068，即表示编号为1，2，…，12的单料烟叶在叶组中所占的百分比例分别为1.43%，12.14%，9.52%，6.7%，23.08%，8.11%，10.1%，3.11%，1.65%，3.7%，19.52%，0.68%。然后，配方设计师就可按照每种单料烟叶的百分比例制作叶组配方来进行抽吸评定，如果抽吸评定的结果不符要求，则将 , , , , 重新取值再进行计算，直到找到满意且符合生产要求的叶组配方为止。

第4章二次型【二次型简介】非线性问题广泛存在于各个科学技术领域，而某些非线性问题在一定的条件下可以转化为线性问题来进行研究。方法之一是通过矩阵的方法将二次型化为标准形，具体包括合同变化法和正交变换法。

第4章二次型验证性实验二次型及标准形【实验目的】掌握二次型及其矩阵表示了解二次型秩、二次型的标准形的概念会用正交变换等方法化二次型为标准形理解正定二次型、正定矩阵的概念，并掌握其判别法【实验要求】掌握分数数据格式format rat、计算矩阵的秩rank、求矩阵的特征值和特征向量eig、单位阵eye等命令

第4章二次型【实验内容】 1．求二次型的矩阵和二次型的秩。 2．用合同变换将二次型化为标准形。 3．用正交变换将二次型化为标准形。 4．判断下列二次型是否正定或负定。（1）（2）

第4章二次型【实验过程】 1．>> format rat >>A=[1 -3/2 -1;-3/2 2 1;-1 1 3] 运行结果： A = 1 -3/ / >>rA=rank(A) rA = 3

第4章二次型 2．>> format rat >> A=[ ; ; ; ]; >> E=eye(4); >> AE=[A,E] 运行结果： AE = >> AE(1,:)=AE(1,:)+AE(2,:); >> AE(:,1)=AE(:,1)+AE(:,2)

第4章二次型运行结果： AE = >> AE(2,:)=AE(2,:)-1/2*AE(1,:); >> AE(:,2)=AE(:,2)-1/2*AE(:,1) 运行结果: / /2 1/

第4章二次型 >> AE(3,:)=AE(3,:)-2*AE(2,:); >> AE(:,3)=AE(:,3)-2*AE(:,2); >> AE(4,:)=AE(4,:)+2*AE(2,:); >> AE(:,4)=AE(:,4)+2*AE(:,2) 运行结果： AE = / /2 1/

第4章二次型 >> AE(4,:)=AE(4,:)+1/2*AE(3,:); >> AE(:,4)=AE(:,4)+1/2*AE(:,3) 运行结果： AE = / /2 1/ /2 -1/2 1/2 1/2 1 得，即经过可逆线性变换，原二次型化为标准形

第4章二次型 3．>> A=[ ; ; ; ]; >> [kesai,lamda]=eig(A) 运行结果： kesai = lamda =

第4章二次型得即正交变换将原二次型化为标准形

第4章二次型 4．（1）>> A=[-2 1 1;1 -6 0;1 0 -4]; >>a1=det(A(1,1)) 运行结果： a1 = -2 >>a2=det(A(1:2,1:2)) a2 = 11 >>a3=det(A(1:3,1:3)) a3 = -38 由上可知，A的奇数阶顺序主子式小于0，而偶数阶顺序主子式大于0，故此二次型是负定的。

第4章二次型（2）>> A=[ ; ; ; ]; >>a1=det(A(1,1)) 运行结果： a1 = 1 >>a2=det(A(1:2,1:2)) a2 = 2

第4章二次型 >>a3=det(A(1:3,1:3)) 运行结果： a3 = 6 >>a4=det(A(1:4,1:4)) a4 = 24 由上可知，A的各阶顺序主子式都大于0，故此二次型是正定的。

第4章二次型设计性实验实验一三元二次方程的三维图形判定【实验目的】 1．从实际问题出发，建立对应的二次型
2．掌握二次型的特性和矩阵 3．掌握空间解析几何中关于二次方程图形的形式【实验要求】掌握Matlab中关于二次型转化的各种命令【实验内容】现给定一个三元二次方程，形如：

第4章二次型如何判定其在三维空间的图形？【实验方案】利用矩阵的特征值和特征向量求解。【实验过程】
把中的变量组成部分转换为二次型，得到：接下来我们需要把它转化成标准形，再进行判定。不妨假设问题为 =C的三维空间图像判定

第4章二次型首先先写出二次型的实对称矩阵在Matlab编辑器中建立M文件如下： A=[2 -2 0;-2 1 -2;0 -2 0];
[P,D]=schur(A) syms y1 y2 y3 y=[y1;y2;y3]； f=[y1 y2 y3]*D*y

第4章二次型运行结果： f =2*y1^2-y2^2+4*y3^2 即通过空间解析几何的知识我们很容易判定：

第4章二次型实验二线性规划方法建模【实验目的】 1．掌握线性规划模型 2．掌握线性规划解的基本理论 3．掌握线性规划的求解方法
【实验要求】学会使用MTLAB软件的linprog命令求解线性规划模型【实验内容】某车间有两台机床甲和乙，可用于加工三种工件。假定这两台机床的可用台时数分别为700和800，三种工件的数量分别为300、500和400，且已知用三种不同机床加工单位数量的不同工件所需的台时数和加工费用（如表4-1所示），问怎样分配机床的加工任务，才能既满足加工工件的要求，又使总加工费用最低？

第4章二次型线性规划模型设产品产量为，称之为决策变量，所得的利润为z，则要解决的问题的目标是使得（总利润）函数有最大值．决策变量所受的约束条件为

第4章二次型问题可归结为求目标函数在约束条件下的最大值问题．目标函数和约束条件都是决策变量的线性函数，即有下面的线性规划模型目标函数：
约束条件：一般地，如果问题的目标函数和约束条件关于决策变量都是线性的，则称该问题为线性规划问题，其模型称为线性规划模型．我们规定线性规划模型的标准型为

第4章二次型对于非标准型的线性规划模型都可以化为标准型，其方法如下：（1）目标函数为最小化问题：令，则；
（1）目标函数为最小化问题：令，则；（2）约束条件为不等式：对于不等号“ ”的约束条件，则可在“ ”的左端加上（或减去）一个非负变量（称为松弛变量）使其变为等式．（3）对于无约束的决策变量：譬如，则令，使得，代入模型即可．线性规划的求解方法 MATLAB中，线性规划问题（Linear Programming）的求解使用的是函数linprog。

第4章二次型表4-2 函数linprog的使用格式

第4章二次型说明：若exitflag>0表示函数收敛于解x，exitflag=0表示超过函数估值或迭代的最大数字，exitflag<0表示函数不收敛于解x；若lambda=lower 表示下界lb，lambda=upper表示上界ub，lambda=ineqlin表示不等式约束，lambda=eqlin表示等式约束，lambda中的非0元素表示对应的约束是有效约束；output=iterations表示迭代次数，output=algorithm表示使用的运算规则，output=cgiterations表示PCG迭代次数。

第4章二次型【实验过程】设在甲机床上加工工件1、2和3的数量分别为x1、x2和x3，在乙机床上加工工件1、2和3的数量分别为x4、x5和x6。根据三种工种的数量限制，有 x1+x4=300 (对工件1) x2+x5=500 (对工件2) x3+x6=400 (对工件3) 再根据机床甲和乙的可用总台时限制，可以得到其它约束条件。以总加工费用最少为目标函数，组合约束条件。

第4章二次型首先输入下列系数： >>f = [13;9;10;11;12;8];
>>A = [ ]; >>b = [700; 800]; >>Aeq=[ ]; >>beq=[ ]; >>lb = zeros(6,1); >>[x,fval,exitflag,output,lambda] = linprog(f,A,b,Aeq,beq,lb);

第4章二次型运行结果： Ans=x = 0.0000 Ans=fval = e+004 Ans=exitflag =1 可见，在甲机床上加工500个工件2，在乙机床上加工300个工件1、加工400个工件3可在满足条件的情况下使总加工费最小。最小费用为11000元。收敛正常。

第4章二次型实验三储层孔隙度的准确计算【实验目的】 1．掌握最小二乘法的原理的基础 2．根据实际工程问题建立适合的数学模型，
【实验要求】掌握MATLAB工箱函数中的最小二乘法函数求解。【实验内容】储层孔隙度的准确计算是油层识别、含油饱和度计算和精细评价的基础。计算地层孔隙度的方法很多，包括声波测井响应方程、交会图、利用岩心分析与声波时差建立回归关系式等方法。但由于储层非均质性强，地层中子孔隙度及密度测井受地层水矿化度、原油密度及岩性的影响比较大，单纯利用声波时差计算孔隙度忽略了密度和中子测井对孔隙度的影响。因此综合采用多测井曲线,利用统计方法建立测井曲线与孔隙度的模型,往往会取得较好的应用效果。

第4章二次型本实验用测井曲线:密度（DEN）、声波（AC）、中子（CNL）、自然伽玛（GR）与岩心孔隙度建立模型，再用最小二乘法求解。数据如表4-4：

第4章二次型【实验方案】 (一)最小二乘法的原理
求解线形方程组时,通常要求未知数的个数与方程式的个数相等,如果方程式的个数多于未知数的个数,就往往无解,这样的方程组称为矛盾方程组"最小二乘法是用来解矛盾方程组的一个常用方法。设有矛盾方程组 (1)

第4章二次型因为，所以找不到能同时满足这n个方程的解。因此转而寻求在某种意义下的近似解，这种近似解当然不是指对精确解的近似(因为精确并不存在)，而是指寻求各未知数的一组取值，使方程组中各式近似地相等，这就是用最小二乘法解矛盾方程组的基本设想。根据这一设想，构造目标函数并极小化 (2)

第4章二次型称为最小二乘问题。在模型（1）中，当是线性函数时，称为（2）为线性最小二乘问题，当是非线性函数时，称（2）为非线性最小二乘问题。（１）线性最小二乘问题在式（2）中，假设我们可以用矩阵乘积形式来表达线性最小二乘问题

第4章二次型则（２）非线性最小二乘问题在式（2）中若是非线性函数，且存在连续的偏导数。解决这类问题的基本思想：通过一系列线性最小二乘问题求非线性最小二乘问题的解。

第4章二次型【实验过程】 MATLAB的程序清单如下： function [B,U,Q,R,F,SS]=dd1(X,Y,m,n)
for i=1:m M(i)=mean(X(:,i)); X1(:,i)=X(:,i)-M(i); end y=mean(Y); Y1=Y-y; for j=1:m S(i,j)=sum(X1(:,i).*X1(:,j)); Sy(i)=sum(X1(:,i).*Y1);

第4章二次型 Sy=Sy'; M=M'; B=inv(S)*Sy; y1=sum(M.*B); b0=y-y1; B=[b0;B]
for i=1:n YY(i)=sum(B(2:m+1).*X(i,:)')+B(1); Y3(i)=Y(i)-YY(i); end Y2=YY-y; U=sum(Y2.^2) Q=sum(Y3.^2) t=Q/(U+Q); R=sqrtm(1-t) f1=U/m; f2=Q/(n-m-1); F=f1/f2 SS=sqrtm(f2)

第4章二次型运行结果： >>load lsqdata >> dd1(x,y,4,28) B = 29.1414
0.1883 U = Q = R = 0.7506 F = 7.4184 SS = 1.5011

第4章二次型根据运行结果，建立方程如下：说明：(1)lsqdata为数据文件，存放原始数据包含两个矩阵，其中，Y=岩心孔隙度。
（2）分别为数据矩阵，变量数，样本容量，方程系数，检验值。

第4章二次型实验四隧道围岩监测位移前推【实验目的】 1.掌握曲线拟合的方法。
2．能够在Matlab中运用一元线性回归方法（多项式拟合）解决某些实际问题。【实验要求】 1．曲线拟合的方法，一元线性回归（多项式拟合）的算法及实施步骤。 2．Matlab中一元线性回归（多项式拟合）命令polyfit的用法。【实验内容】在隧道开挖过程中，通常要对围岩（隧道周围等于其横剖面中最大尺寸的3—5倍范围之内的岩体）进行位移量的监测，并以此来对围岩以及整条隧道的稳定性和安全性做出评价和预测。但是，在通常情况下，监测点的埋设受各种

第4章二次型因素的影响，并不能完全紧跟掌子面（开挖隧道时对隧道前方岩体进行挖掘的工作面）掘进而进行，围岩在采集监测初值之前已产生部分位移，也就是说监测到的围岩位移量丢失了隧道开挖前期部分重要的围岩位移释放量。因此，这样得到的位移值不是围岩的实际位移值，这会对围岩和隧道稳定性的评价和预测产生很大的影响，甚至有可能导致对围岩和隧道的稳定状态的判定发生偏差，从而引发工程事故。这里就要求我们根据已有的围岩监测位移量估计隧道开挖前期丢失的那部分位移量，对监测到的位移量做出修正以获得真实的围岩位移释放量。某隧道围岩位移监测点的埋设滞后其开挖时间点7天，其中一个监测点前10天的累计位移值如下表所示，要求按照表4-5中数据建立相应的数学模型，对该监测点的位移前期释放量做出估计，并对该监测点的位移值进行修正最终获得围岩的真实位移量。

第4章二次型【实验方案】为了根据已有的位移值估计前面释放的位移值，首先需要做出图形来更加直观地观察围岩的累计值的发展变化趋势，了解位移时间序列的趋势。那么可以以监测时间为横坐标，累计位移值为纵坐标，用MATLAB中的plot函数画出的二维散点图。

第4章二次型若从绘制的监测位移散点分布图中可以看出散点的分布近似于一条的直线，那么我们就可以用一元线性回归的方法来拟和这些散点，求出回归方程，即求出该直线的具体函数表达式，这样的话我们就可以计算出当该直线的自变量取负值时的函数值，即求得监测开始前的位移值，进而对监测到的位移值作出修正。【实验过程】 >> X=[ ]; >> Y=[ ]; >> plot(X,Y,'.'); >> axis([0, 10, 0, 0.8]); >> xlabel('时间（天）'); >> ylabel('累计位移值（mm）');

第4章二次型运行结果：图4-1 监测位移散点图

第4章二次型 >> X=[ ]; >> Y=[ ]; >> P=polyfit(X,Y,1) %对X,Y进行一元线性回归，即多项式拟和运行结果： P = 因此所求得的直线的函数表达式为y=0.0687x 。

第4章二次型由于在该问题中，监测点的埋设滞后其开挖时间点7天，分别取x=-7，-6，-5，-4，-3，-2，-1，0，得到y= ，，，，，，，。然后将当x=-7时的y值置零,即将每个y值和本来已有的监测位移值都加上0.4976，这样就得到了从隧道开挖时刻起围岩的实际位移释放值（取2位小数），如表4-7所示。

第5章综合性实验综合性试验1 原因和可识别性一问题背景与实验目的
求解线性方程组是线性代数中最基本问题，线性代数主要关注的是存在性、唯一性和稳定性问题，而对反问题所涉及的不稳定性讨论得并不多。线性代数的正问题包含确定线性变换的表示矩阵；给定的矩阵和一个n维向量 ,确定m维向量。而寻找所有满足的解的原因反问题可能比其他基本的线性代数问题得到更多的关注。一个很少被提到反问题是识别问题；确定矩阵，使得给定的“输入-输出”对的集合满足。这个模型是实矩阵的反问题的基本表达。

第5章综合性实验原因反问题求向量满足 ,其中是一个给定的实矩阵，。
原因反问题求向量满足 ,其中是一个给定的实矩阵，。存在性问题这个反问题的解存在的充分必要条件是。由线性代数知识，子空间就是所有的列向量的线性组合形成的的子空间，并且确定是否属于（即解是否存在），以及寻找所有解是通过一个很好的算法—高斯消去法来实现的。唯一性问题通过的零空间来解决的。同时，高斯消元法是一个刻画零空间的有效手段，因此可以用来解决唯一性问题。稳定性的解关于右端项的扰动的稳定性，可以用矩阵的条件数来量化。现在，我们分析在怎样一个相对误差内会导致解的相对大的改变。

第5章综合性实验（i）对于每个方程存在唯一解的情形假设是右端项的一个扰动，则相对于的大小扰动的大小可以
用来衡量,其中的范数。令是方程组对应的唯一解，而是方程组的解。那么所以矩阵范数给出了解在右端扰动下所产生的误差的界，且有如下的相对误差估计，因此，其中称为矩阵的条件数。因此条件数给出了右端项相对误差的上界。对

第5章综合性实验于大条件数的矩阵，即病态矩阵，右端相对小的扰动会引起解的相对大的变化。在这种意义下，病态线性方程组是不稳定的。
（ii）对于无解的或有无穷多解的线性方程组当不属于的矩阵的值域时，方程组无解。我们希望引入它的一类广义解，并且这一广义解总是存在的。考虑广义解是一种最小二乘解，即解向量使得在所有的中范数达到最小。如果是的最小二乘解（广义解），那么对于任意向量，函数在时达到最小值，其中为熟悉的欧几里德内积。

第5章综合性实验由最小值的必要条件得到，所以对所有，。因此，如果是最小二乘解，那么，其中是的转置矩阵。
得到，所以对所有，。因此，如果是最小二乘解，那么，其中是的转置矩阵。相反的，如果，那么对于任何，即为的最小二乘解。

第5章综合性实验所以的最小二乘解就是的解。由线性代数理论知，，从而总是有一个解。因此对于任意。所以任意线性方程组有一个最小二乘解。另外，如果是一个最小二乘解，那么对于任意的也是最小二乘解，即最小二乘解的集合构成平行于零空间的超平面。因此，如果有一个非平凡的零空间，即，那么有无限多个最小二乘解。然而其中有一个最小二乘解与零空间正交的，且这样的最小二乘解至多只有一个，因为如果和都是最小二乘解，并且与正交，那么也与正交。而，所以。因此，即。另一方面，这种与零空间正交的最小二乘解总是存在的，所以任意线性方程组总

第5章综合性实验存在唯一的与系数矩阵的零空间正交的最小二乘解。这样，如果我们接受这一广义解的概念，那么每个线性系统都有唯一（广义）解。最后,我们简要的考虑一下识别的问题,即给定关于的向量对，确定矩阵的反问题。对于每个向量对，我们称为输入，而为相应的输出。我们的工作是通过“询问”适当的，观察输出的来识别“黑箱” 。因为我们控制输入，可以安排它们成为线性无关的，用一个矩阵来表示输入，其中矩阵的每一列分别为输入向量。类似的,相应的输出向量可用矩阵表示。如果存在唯一的矩阵满足，那么称可以通过矩阵对来识别。

第5章综合性实验如果，此时是可逆的，则，从而是可识别的。这种情况下，一个经过简单修改的高斯消元法可以成为识别的一个算法。这一方法依赖于这一事实，即如果，那么。因此可以通过高斯-约当算法求解模型矩阵的转置：如果，识别的输入输出信息是不充分的。因为在这种情况下，存在一个向量与正交。令是矩阵（它的第一行是，其他行都是零向量）。那么是的零矩阵。因此，如果，同样的，所以从信息是无法识别的。

第5章综合性实验二实验内容确定矩阵满足，其中 , 。三实验方案等式两边同乘以矩阵的逆，，即
确定矩阵满足，其中 , 。三实验方案等式两边同乘以矩阵的逆，，即在Matlab编辑器中建立M文件如下： X = [1 1 1;-1 1 0;0 1 1]; B = [-2 6 3;-1 0 0]; A = B*inv(X) 运行结果为： A =

第5章综合性实验综合性试验2 断层成像一问题背景与实验目的
断层成像技术的历史背景可以追溯到柏拉图从投到墙上的影子重构实物的问题，一直到柯马克的射线断层摄影术的发展。所有这些问题的本质都是从某个投影中重构实体，这样的问题一般说来是欠定的。也就是说，对于相对贫乏的数据，无法对物体进行唯一重构。然而，重构相应数据的近似解的方法在很多技术中是很有用的。下面我们将讨论简单的集合推动算法，“代数重构技术”（ART），它可以用来求射线断层摄影术中非常基本的线性代数问题的近似解。

第5章综合性实验将ART算法应用到一个很简单的射线断层摄影术问题中。考虑物体为排列的“像素”。我们假定物体在每个像素内是一样的，但像素与像素之间可能存在差异，我们的目标是得到关于像素之间的变化情况的图片。假设物体被通过像素中心的放射线所照射。我们从第一行开始从左到右排列这些像素。光线通过物体的路线定义为观察向量，即一个维的行向量，它的元素取决于光线是否通过该像素而为0（否）或1（是）。以一个为例，在图5-1中光线的观察向量为

第5章综合性实验图5-1 物体的切片假设材料的每个像素吸收射入它的放射线的一部分，从而所吸收的部分是区间上的一个数。如果该数值接近于0说明像素是接近透明的，而接近1的值表明像素基本不透明，我们称第个像吸收部分为该像素的吸收系数。 1 2 3 4 5 6 7 8 9

第5章综合性实验假设一束光线的观察向量为，其中，以及穿过物体的吸收系数为。如果光线穿过第个像素，那么放射线的部分被吸
其中，以及穿过物体的吸收系数为。如果光线穿过第个像素，那么放射线的部分被吸收，而部分传递出去了。从物体中穿透出去放射线的部分为其中为光线的观察向量，或等价的，注意到每个为负数，所以我们可以把这一关系写成

第5章综合性实验其中为正数。因此这一射线断层摄影技术问题等价于寻找由个观察向量和相关的右端项确定的个超平面的交集中的一个点。此外，对于某一特定观察，其右端项是射线穿过的像素相应的之和。给定个观察向量和相应的值，我们把这些观察值作为观察矩阵的行向量，其第个元素为，并且把与这些观测值相应的测量值作为维行向量。给定维行向量表示初始近似值，对于任何自然数，执行了次ART 算法的投影循环后可以得到后量的一个近似值。由程序“art1”（调用格式 :）是用特别的扭曲来实现，即每次循环后，将近似解投影到

第5章综合性实验平面（对于所有的）由“art1”得到的近似解可以通过程序“displa”表示出来，其结果为的网点，每个网点产生适当的阴影。下面是简单的摄影断层扫描解码技术的例子。考虑一个的物体，有一个几乎全黑的主对角线和一个灰色的交错对角线。我们可以将它编码，如图5-2所示，用向量表示这个物体图5-2 样品 1 0.5

第5章综合性实验假设我们有物体的六个观察值，包括三个行片层，沿主对角线向下，沿另一根对角线向上和沿中间列向下。那么观察矩阵为
且相应的测量向量为

第5章综合性实验二实验内容一条虫子在一个划分为网格的苹果上打洞。穿过行（从上到下）的水平探测测量值分别为0，2，1，2，0，而沿列（从左到右）垂直探测测量值分别为0，1，3，1，0.主对角线（从西北到东南）之和为1.而另一条对角线（西南到东北）之和为3.使用零向量为初始近似，调用程序“art1”和“displa”计算虫子的图像。三实验内容在Matlab中分别建立以下四个M文件

第5章综合性实验 %求解迭代后的逼近值 function [x]=art1(V,mu,x0,n); m=size(V,1);
xold=x0; for j=1:n for k=1:m; x=proj(V(k,:),xold,mu(k)); xold=max(x,zeros(1,size(x,2))); x=xold; end;

第5章综合性实验 %求解向量在某一平面上的投影 function [p]=proj(v,x,mu); normv=norm(v);
if normv==0 normv=1; end; p=x+((mu-x*v')/normv^2)*v;

第5章综合性实验 %生成网格点并对其进行填充 function [n]=displa(x) axis('equal');
n=sqrt(size(x,2)); xmax=max(x); for k=1:n; for l=1:n; [X,Y,s]=shade(l-1,n-k,x(n*(k-1)+l)/xmax); fill(X,Y,s); hold on; end;

第5章综合性实验 %填充效果 function [X,Y,s]=shade(i,j,p); X=i+[0 1 1 0];
Y=j+[ ]; s=(1-p)*[1 1 1];

第5章综合性实验综合性试验3 艾滋病疗法的评价及疗效的预测【问题的提出】
艾滋病是当前人类社会最严重的瘟疫之一，医学全名为“获得性免疫缺损综合症”，英文简称AIDS，它是由艾滋病毒（医学全名为“人体免疫缺损病毒”, 英文简称HIV）引起的。这种病毒破坏人的免疫系统，使人体丧失抵抗各种疾病的能力，从而严重危害人的生命。人类免疫系统的CD4细胞在抵御HIV的入侵中起着重要作用，当CD4被HIV感染而裂解时，其数量会急剧减少，HIV将迅速增加，导致AIDS发作。艾滋病治疗的目的，是尽量减少人体内HIV的数量，同时产生更多的CD4，至少要有效地降低CD4减少的速度，以提高人体免疫能力。

第5章综合性实验迄今为止人类还没有找到能根治AIDS的疗法，只是通过一些药物疗法来延缓AIDS的发病时间，但是由于药物的副作用很大和费用较高，所以应该选择适合的疗法来治疗。现在有美国艾滋病治疗实验机构ACGR公布的两组药物治疗数据同时服用zidovudine（齐多夫定），lamivudine（拉美夫定）和indinavir（茚地那韦）3种药物的300多名病人每隔几周测试的CD4和HIV的浓度（每毫升血液里的数量）。和将1300多名病人随机地分为4组，每组按下述4种疗法中的一种服药，大约每隔8周测试的CD4浓度（这组数据缺HIV浓度，它的测试成本很高）。4种疗法的日用药分别为：600mg zidovudine或400mg didanosine（去羟基苷），这两种药按月轮换使用；600 mg zidovudine加2.25 mg zalcitabine（扎西他滨）；600 mg zidovudine加400 mg didanosine；600 mg zidovudine加400 mg didanosine，再加400 mg nevirapine（奈韦拉平）。

第5章综合性实验我们通过对以上问题的分析完成了以下问题的解：
（1）利用数据，预测继续治疗的效果，确定最佳治疗终止时间（继续治疗指在测试终止后继续服药，如果认为继续服药效果不好，则可选择提前终止治疗）。（2）利用数据，以CD4为标准评价了四种疗法的优劣，并对较优的疗法预测继续治疗的效果，并确定最佳治疗终止时间。（3）考虑艾滋病的价格因素，对于不发达国家的病人来说，（2）中的预测和评价又有了相应的改变。

第5章综合性实验【问题的分析】 1、对问题一的分析：
问题要求利用附件1的数据，预测继续治疗的效果，或者确定最佳治疗终止时间。数据ACTG320是使用同时服用3种药物的治疗方法得到的300多名病人每隔几周测试的CD4和HIV的浓度。把附件1的扩展名txt改成xlt（即电子表），运用excel可统计出结果为355位病人的1665组检测数据，其中病人、时间、CD4数量和HIV的浓度是要考虑的四个量。以下是处理数据的原则。 ⑴ 剔除缺少HIV的数据运用excel统计，有100组数据缺少HIV数值，仅占总量的6%，为了便于研究，首先把这些组数据剔除掉。

第5章综合性实验 ⑵ 利用差值法把数据组时间统一到第0、4、8、24、40周
将1665组检测值对时间作散点图（见图1），从图1可以发现：图形大致呈折线形变化，且在第0、4、8、24、40周附近呈水平，为了简化对问题的处理，可在总体上只考虑这五周CD4的大小和HIV的数值（以下简记为CD4和HIV）来进行总体拟合，对恰巧不全是这几周检测的，我们可按下面四条原则行处理： ①按相邻原则由左右两测点推测中间点。我们假设病人的病情在两个测试中间段的变化是平稳的，具体办法：如已知第0、3、8、23、40周，可由第3周与第8周的CD4和HIV分别按直线求出第4周的CD4和HIV，同理可由第23周与第40周的CD4和HIV分别按直线求出第24周的CD4和HIV。

第5章综合性实验 ②按就近原则由已知最近两点推测端点。具体办法：如已知第1、4、8、24、38周，可由第1周与第4周的CD4和HIV分别按直线求出第0周的CD4和HIV，同理可由第24周与第38周的CD4和HIV分别按直线求出第40周的CD4和HIV。 ③按多余点舍弃原则。具体办法：如已知第0、4、8、24、41、50周，可由第24周与第41周的CD4和HIV分别按直线求出第40周的CD4和HIV，对多余的第50周可舍弃。 ④按缺少点不补原则。具体办法：如已知第0、4、8、24周，对缺少第40周CD4和HIV，不再补充，在求总体均值时要注意到总量的多少。 ⑶第0周情况可作为病情分类的依据。由于第0周是病人治疗前的检测结果，这既可作为病情分类的依据，也可作为治疗结果的参照。当然，在对总体拟合是也是不可缺少的点。

第5章综合性实验 2、对问题二的分析：问题要求利用附件2的数据，仅以CD4为标准来评价4种疗法的优劣，并对较优的疗法预测继续治疗的效果，或者确定最佳治疗终止时间。附件2涉及到五个量：病人、疗法、年龄、时间、Log(CD4数+1)。评价4种疗法的优劣时，首先按不同的治疗法将数据分成四组，然后对每组的CD4分别按问题一的方法进行处理，统计出各组CD4值增加的个数，据此运用统计概率模型可以求出四种治疗法的有效治疗率，由有效治疗率大小即可以判别四种疗法的优劣。对较优的疗法，然后仿照问题一依据病人治疗前的CD4值将病人按病情分成早期、中期和晚期三类，再依据年龄将病人分成青壮年和中老年，这样交叉可得六类病人，最后仿照问题一求出各类数据拟合函数曲线图，据此可以预测继续治疗的效果，及确定最佳治疗终止时间。

第5章综合性实验 3、对问题三的分析：近年来，由于不发达国家因为种种原因，感染艾滋病人数很多，而经济又相对落后，治疗费用是必须考虑的一大重要因素，为此，对问题二的模型的要进行改进，要把几种药品价格加入模型进行分析与处理。对个体而言，由于不同的疗法的费用是固定的，我们要从总体上进行考虑，在疗效相同的条件下，对每种疗法进行实验设计，然后可求得每种费用在相同时间内的花费多少，从而可以判四种疗法的优劣。对最佳停药时间，应该为有费用限制出现在无费用限制之后。

第5章综合性实验【模型的假设】 1、假设病人的病情在两个测试周中间时段的变化是平稳的；
2、为了方便起见，规定早期、中期、晚期编排序号为1、2、3； 3、规定第0、4、8、24、40周编排序号为第0、1、2、3、4次； 4、每月按28天计算即每个月4周； 5、对仅有初始检测值的视为未进行治疗，在研究时可剔除； 6、所有数据均为原始数据，来源真实可靠。

第5章综合性实验【模型的建立与求解】 1、模型Ⅰ 数据拟合模型
本模型从CD4数量越多则病人的免疫力越强，HIV的数值越小则病人的病情越轻出发，建立了健康指标模型，给出每个病人每个时段的健康指标，依据病人治疗前的健康指标将病人按病情分成早期、中期和晚期三类，然后利用最小二乘法原理分别求出各类数据拟合函数曲线图，据此可以预测继续治疗的效果，及确定最佳治疗终止时间。 2、模型Ⅱ 统计概率模型模型本模型首先对四类治疗法分成组，然后对每组的CD4分别按问题一的方法进行处理，统计出各组CD4值增加的个数，据此运用统计概率模型可以求出四种治疗法的有效治疗率，由此可以判别四种疗法优劣。

第5章综合性实验 3、模型Ⅲ 数据拟合模型本模型针对模型Ⅱ求出的较优的疗法，然后仿照模型Ⅰ依据病人治疗前的CD4值将病人按病情分成早期、中期和晚期三类，再依据年龄将病人分成青壮年和中老年，这样交叉可得六类病人，最后仿照模型Ⅰ求出各类数据拟合函数曲线图，据此可以预测继续治疗的效果，及确定最佳治疗终止时间。 4、模型Ⅳ 统计实验模型对个体而言，由于不同的疗法的费用是固定的，我们要从总体上进行考虑，为此，本模型是从总体上在疗效相同的条件下，对每种疗法设计进行实验设计，然后可求得每种费用在相同时间内的花费多少，从而可以判四种疗法的优劣。对最佳停药时间，应该为有费用限制出现在无费用限制之后。

第5章综合性实验一、问题一的分析与求解 1、对问题的分析由于利用附件1的数据，预测继续治疗的效果，或者确定最佳治疗终止时间。
⑴相关量的定义 ①健康指标为了能更好地反映出同时服用zidovudine、lamivudine和indinavir这三种药物的治疗效果，由问题的分析我们知道，CD4数量越多，病人的免疫力越强，HIV的数值越小则病人的病情越轻，为此兼顾这两个方面，下面我们定义一个新的指标函数。定义1 在同一时刻，病人血液中的CD4浓度与HIV浓度的比值反映一个人的健康状况，我们称为健康指标，记作k。

第5章综合性实验由定义，我们可以给出第i个病人在第j次检测时的健康指标为
其中、分别表示第i个病人在第j次检测时的CD4与HIV浓度。因为在研究问题过程中单位对结果无影响，我们对健康指标的单位可以不予考虑。显然，一个病人的健康指标值越大表示这个人的病情越轻，反之病情越重。 ②病情好转速度为了确定最佳治疗终止时间，下面我们由健康指标与时间二者的相对改变量给出以下定义。定义2 我们称单位时间（周）内一个病人健康指标的变化量为其病情好转速度，记作v。

第5章综合性实验由定义，我们可以给出第i个病人在第j次检测时的病情好转速度为
其中表示第i个病人在第j次检测时的健康指标，表示第i个病人进行第j次检测的周数。若为正值，表示艾滋病人在服用这三种药物后病况好转，反之，若为负值，表示艾滋病人在服用药物后病情趋于恶劣，即此时再服用这三种药物，对病人将产生负作用。综合考虑每个病人各时期的病情好转速度，我们判断每个病人是继续服药还是提前终止服药并确定最佳治疗终止时间。

第5章综合性实验 ⑴模型的准备 ①对病期的划分
2、模型Ⅰ 数据拟合模型 ⑴模型的准备 ①对病期的划分由问题的分析，可由初始检测结果对病情进行分类，CD4 绝对数≥500个/ vl时为无症状HIV感染期，我们记为早期； CD4 绝对数为 /vl时为有症状感染期，我们记为中期；CD4 绝对数＜200/vl时为AIDS 期，我们记为晚期。考虑到附件1中CD4浓度已乘以0.2个/vl，为保持一致，将上述CD4的数目乘以0.2，可以得出与附件1中的数据单位一致的艾滋病病期划分标准（见表5-3）。表5-3 艾滋病病期划分标准（一）病期晚期中期早期修正CD4值范围 0~40 40~100 100~

第5章综合性实验 ②各病期人数的分布利用EXCEL软件对数据进行自动筛选，得到最初开始测试CD4即还未开始服用这三种药物时336名患者在各病期的分布情况（见表5-4）。表5-4 治疗前艾滋病人在各病期所占的比例患病程度频数m比例r晚期 %中期 %早期 %③统一求出第0、4、8、24、40周的CD4和HIV浓度我们利用问题分析中的具体处理办法，利用差值法把时间统一到第0、4、8、24、40周，求出每个病人在第0、4、8、24、40周的测试的CD4和HIV浓度。患病程度频数m 比例r 晚期 113 33.63% 中期 94 27.98% 早期 129 38.39%

第5章综合性实验 ④各期平均健康指标要预测病人继续治疗的效果，或者确定最佳治疗终止时间，总体健康指标是解决解决问题的关键，在统计学上总体健康指标往往用样本的平均指标来估计，因此下面我们先求平均健康指标。结合表5-4中三个病期的病人数，可建立以下平均健康指标模型，其中分别表示早、中、晚期，表示t期第i个病人在第j次的健康指标。利用EXCEL软件从附录3的数据中分别筛选出个各病期病人的五个时段的健康指标值，然后对上述模型用Matlab软件编程求解（见附录4）。求得的各个病期病人五个检测时间的平均健康指标结果见表5-5。 ⑵模型的建立与求解 1）数据的拟合原理：最小二乘法原理

第5章综合性实验表5-5 各期各检测时间病人的平均健康指标根据表3中已有的各病期的5个检测时刻的健康指标数据，为了方便起见，取为，其中依次为0、4、8、24、40五个点。在以平面直角坐标系中我们很容易作出各自的散点图，通过观察各散点图中散点的变化趋势，我们可构造含有一个或若干个未知参数的拟合函数，拟合函数可以为直线、抛物线、指数函数或其它函数，这取决于经验或软件。对拟合的优劣，可根据最小二乘法原则，拟合函数使得函数在点处的函数值与观测数据偏差的平方和达到最小。检测时间病期 4 8 24 40 早期中期晚期

第5章综合性实验二、问题二的分析与求解 1、对问题的分析
问题要求利用附件2的数据，仅以CD4为标准来评价4种疗法的优劣，并对较优的疗法预测继续治疗的效果，或者确定最佳治疗终止时间。附件2涉及到五个量：病人、疗法、年龄、时间、Log(CD4数+1)。评价4种疗法的优劣时，首先我们按疗法把病人分开成四组，然后对每组的CD4分别按问题一的方法进行处理，再对各种疗法进行排序；其次我们对每组再按年龄分成若干组进行讨论，再细考虑各种疗法对不同年龄组的疗效的好坏，以期得到较优的结果。

第5章综合性实验 2、模型Ⅱ 统计概率模型模型 ⑴模型的分析
根据对文献资料的查阅，要判断一种疗法治疗效果的好坏，主要是通过以下三个方面进行评估：病毒学指标（即HIV浓度指标）、免疫学指标（即CD4浓度指标）和临床症状，其中病毒学的改变是最重要的指标。但是因为HIV浓度测试的成本很高，而且在本题中未给出临床症状的指标，所以在模型中我们仅从免疫学指标即CD4浓度方面来评价这四种疗法的优劣。 ⑵模型的准备首先我们利用EXCEL统计软件从附件2给出的数据中分别将按照疗法1、疗法2、疗法3、和疗法4的病人的CD4与对应的测试时刻时间整理出来（各疗法对应的数据见附录EXCEL表格“问题2”），为了评价4种疗法的优劣，我们以

第5章综合性实验各疗法的治疗效果为标准，并从每个病人最后一次测试时的CD4浓度2与第一次测试时的CD4浓度的变化考虑治疗效果，用表示第i种疗法的治疗效果：其中与分别表示用第i种疗法第一次和最后一次测试病人的CD4浓度。若，表示用第i种疗法治疗对病人有效，而且值越大疗效越好；若，则表示此时用第i种疗法治疗艾滋病将对病人产生负作用。 ⑶模型的建立运用电子表EXCEL，我们可以准确地统计出每种疗法中治疗有效与无效的病例数，具体结果见表5-6。

第5章综合性实验令分别表示第i种疗法有效个数和总个数，则对各种疗法的有效治疗率也即总的疗效可建立如下统计概率模型： ⑷模型的结果
表5-6 四种疗法有效无效治疗个数分布表令分别表示第i种疗法有效个数和总个数，则对各种疗法的有效治疗率也即总的疗效可建立如下统计概率模型： ⑷模型的结果利用表5-6中的数据和概率模型，运用MATHEMATICA软件很容易求各种疗法对应的有效治疗率（见表5-7）。治疗方案疗法1 疗法2 疗法3 疗法4 有效个数 88 124 123 331 无效个数 197 168 174 282 总个数 285 292 297 613

第5章综合性实验结论：由有效治疗率的大小可知， =0.54最大，也即疗法4的治疗方案最好。表5-7 各疗法对应的有效治疗率治疗方案
表5-7 各疗法对应的有效治疗率结论：由有效治疗率的大小可知， =0.54最大，也即疗法4的治疗方案最好。治疗方案疗法1 疗法2 疗法3 疗法4 有效治疗率 0.31 0.42 0.41 0.54

第5章综合性实验 3、模型Ⅲ 数据拟合模型 ⑴模型的准备 ①对病期的划分
3、模型Ⅲ 数据拟合模型 ⑴模型的准备 ①对病期的划分问题一中病期的划分为：早期CD4值范围为100以上，中期CD4值范围为40~100，晚期CD4值范围为0~40。由于附录2中的CD4的浓度是经过Log(CD4+1)处理所得的值，通过观察可以发现其范围在0~6.2971。为了便于对病期的划分，我们要对问题一中病期的划分的CD4作相应的处理。由Log41≈3.7136，Log 101≈4.6151，可相应得到具体的病期划分标准（见表5-8）。表5-8 艾滋病病期划分标准（二）病期晚期中期早期 Log(CD4+1)值的范围 0~3.7136 3.7136~4.6151 4.6151~

第5章综合性实验 ②各病期人数的分布按照表5-8给出的艾滋病病期划分标准，对剔除异常值的数据利用EXCEL软件进行自动筛选，得到最初开始测试CD4即还未开始服用这三种药物时336名患者在各病期的分布情况（见表5-9）。表5-9 治疗前艾滋病人在各病期所占的比例患病程度年龄段频数m 比例r 晚期青壮年 161 中老年 85 中期 34 19 早期 5 4

第5章综合性实验 ③统一求出第0、8、16、24、32、40周的CD4浓度
结合表8中三个病期的病人数，可建立以下平均CD4浓度模型，早期两个段在各检测时间的CD4浓度的平均值：其中分别表示早、中、晚期，表示t期第i个病人在第j次的CD4浓度。

第5章综合性实验利用EXCEL软件从附录5的数据中分别筛选出个各病期不同年龄段病人的五个时段健康指标，再由上述模型用Matlab编程求解。各个病期不同年龄段病人五个检测时间的平均CD4浓度，结果见表5-10。表5-10 不同阶段在不同检测时间的CD4的均值 8 16 24 32 40 早期青壮年 3.4393 16.628 30.288 55.728 9.9985 中老年 254 中期晚期

第5章综合性实验三、问题三的分析与求解 1、对问题的分析
在不发达的国家和地区，人们收入普遍低下，大多数艾滋病患者的个人经济能力有限，而接受治疗需要支付高额的费用，这是患者难以承受的。尽管国家有政策对艾滋病患者进行一定帮助，但帮助往往是有限的，不是解决问题的根本。由于患者无法得到及时和有效的治疗，其病情会迅速恶化。现实生活中费用和疗效之间往往是对立的，即疗效高的疗法所花费的费用就高，而费用低的疗效相对要差得多。因此，如何能够花费少又能得到较好的疗效，这对不同经济承受力的人是值得考虑的一个有意义的事。我们在处理这个问题时，重点兼顾到费用和疗效这一对立问题的两个方面。问题三需要在这两个因素中找出一个合理的匹配关系，使得两个方面的都得到较好的兼顾。

第5章综合性实验 2、模型Ⅳ 统计实验模型 ⑴模型的分析
为了解决问题三，下面我们来设计一个统计实验。要对四种方法进行优劣比较，先固定一个量（比如各疗法的有效人数，由问题二记为），由前面的统计概率模型：我们有其中表示第i种疗法的疗效，然后比较个人在一段时间内所花费的费用，显然，费用最小者为较优的治疗方法。这可由下面单目标模型来求解。 ⑵模型的建立取第i种疗法为一个人在一段时间内的全花费，对第i种疗法的个人在一段时间内所花费的费用可由下面公式给出则较优的治疗方法的花费为：

第5章综合性实验 ⑶模型的求解 ①确定样本容量
我们任意选取中一个值，如选，又由表6中四种疗法对应的有效治疗率依次为 0.31、 0.42、 0.41、 0.54。 ②确定治疗期限可根据具体情况进行确定。在这里为了方便疗法1两种药的按月交替使用，故选t=8周（二个月）时间进行比较。由4种疗法的日用药分别为：600mg zidovudine或400mg didanosine（去羟基苷），这两种药按月轮换使用；600 mg zidovudine加2.25 mg zalcitabine（扎西他滨）；600 mg zidovudine加400 mg didanosine；600 mg zidovudine加400 mg didanosine，再加400 mg nevirapine（奈韦拉平）。

第5章综合性实验疗法1的费用为：28×1.60+28×0.85=70.3(美元);
又由600mg zidovudine 1.60美元，400mg didanosine 0.85美元，2.25 mg zalcitabine 1.85美元，400 mg nevirapine 1.20美元。则疗法1的费用为：28× ×0.85=70.3(美元); 疗法2的费用为：56×( )=193.2(美元); 疗法3的费用为：56×( )=137.2(美元); 疗法4的费用为：56×( )=204.4(美元). 则每人各疗法在t=8周内的费用如表5-11：表5-11：各疗法在二个月内的费用疗法疗法1 疗法2 疗法3 疗法4 费用（美元） 70.3 193.2 137.2 204.4

表5-12：第i种疗法的个人在8周内所花费的总费用
第5章综合性实验 ③确定个人在一定时间内的总费用。利用公式得到在二个月内的总费用（见表5-12）。表5-12：第i种疗法的个人在8周内所花费的总费用从表5-12可以看出，费用由少到大依次为：疗法1、疗法3、疗法4、疗法2。综合考虑疗法1效果最好，疗法2效果最差。 ⑷多次实验验证：对疗法1的人数取不同的值，或治疗期间t取不同的时段，同样我们可以得到的相同的结论，即疗法1效果仍然最好，疗法2效果最差。疗法疗法1 疗法2 疗法3 疗法4 人数 1000 738 756 574 总费用 70300 142580 103720 117330

第5章综合性实验 ⑸结论：与问题二比较，疗法优劣改变了，综合考虑疗效、费用，我们得出疗法1的效果好。由于疗法的改变，他们的预测（或者提前结束）必然改变，再加上费用的限制，最佳停药时间一定在无费用限制的情况下出现的早。

第5章综合性实验综合性试验4 电力市场输电阻塞管理模型【问题的提出】
我国电力系统的市场化改革正在积极、稳步地进行。2003年3月国家电力监管委员会成立，2003年6月该委员会发文列出了组建东北区域电力市场和进行华东区域电力市场试点的时间表，标志着电力市场化改革已经进入实质性阶段。可以预计，随着我国用电紧张的缓解，电力市场化将进入新一轮的发展，这给有关产业和研究部门带来了可预期的机遇和挑战。电力从生产到使用的四大环节——发电、输电、配电和用电是瞬间完成的。我国电力市场初期是发电侧电力市场，采取交易与调度一体化的模式。电网公司在组织交易、调度和配送时，必须遵循电网“安全第一”的原则，同时要制订一个电力市场交易规则，按照购电费用最小的经济

第5章综合性实验目标来运作。市场交易-调度中心根据负荷预报和交易规则制订满足电网安全运行的调度计划――各发电机组的出力（发电功率）分配方案；在执行调度计划的过程中，还需实时调度承担AGC（自动发电控制）辅助服务的机组出力，以跟踪电网中实时变化的负荷。设某电网有若干台发电机组和若干条主要线路，每条线路上的有功潮流（输电功率和方向）取决于电网结构和各发电机组的出力。电网每条线路上的有功潮流的绝对值有一安全限值，限值还具有一定的相对安全裕度（即在应急情况下潮流绝对值可以超过限值的百分比的上限）。如果各机组出力分配方案使某条线路上的有功潮流的绝对值超出限值，称为输电阻塞。当发生输电阻塞时，需要研究如何制订既安全又经济的调度计划。

第5章综合性实验电力市场交易规则： 1. 以15分钟为一个时段组织交易，每台机组在当前时段开始时刻前给出下一个时段的报价。各机组将可用出力由低到高分成至多10段报价，每个段的长度称为段容量，每个段容量报一个价（称为段价），段价按段序数单调不减。在最低技术出力以下的报价一般为负值，表示愿意付费维持发电以避免停机带来更大的损失。 2. 在当前时段内，市场交易-调度中心根据下一个时段的负荷预报,每台机组的报价、当前出力和出力改变速率，按段价从低到高选取各机组的段容量或其部分（见下面注释），直到它们之和等于预报的负荷，这时每个机组被选入的段容量或其部分之和形成该时段该机组的出力分配预案（初始交易结果）。最后一个被选入的段价（最高段价）称为该时段的清算价，该时段全部机组的所有出力均按清算价结算。

第5章综合性实验市场交易-调度中心在当前时段内要完成的具体操作过程如下：
1.监控当前时段各机组出力分配方案的执行，调度AGC辅助服务，在此基础上给出各机组的当前出力值。 2.作出下一个时段的负荷需求预报。 3.根据电力市场交易规则得到下一个时段各机组出力分配预案。 4.计算当执行各机组出力分配预案时电网各主要线路上的有功潮流，判断是否会出现输电阻塞。如果不出现，接受各机组出力分配预案；否则，按照如下原则实施阻塞管理：输电阻塞管理原则： (1)调整各机组出力分配方案使得输电阻塞消除。 (2)如果(1)做不到，还可以使用线路的安全裕度输电，以避免拉闸限电（强制减少负荷需求），但要使每条线路

第5章综合性实验上潮流的绝对值超过限值的百分比尽量小。
(3)如果无论怎样分配机组出力都无法使每条线路上的潮流绝对值超过限值的百分比小于相对安全裕度，则必须在用电侧拉闸限电。 (4)当改变根据电力市场交易规则得到的各机组出力分配预案时，一些通过竞价取得发电权的发电容量（称序内容量）不能出力；而一些在竞价中未取得发电权的发电容量（称序外容量）要在低于对应报价的清算价上出力。因此，发电商和网方将产生经济利益冲突。网方应该为因输电阻塞而不能执行初始交易结果付出代价，网方在结算时应该适当地给发电商以经济补偿，由此引起的费用称之为阻塞费用。网方在电网安全运行的保证下应当同时考虑尽量减少阻塞费用。

第5章综合性实验需要做的工作如下： 1.某电网有8台发电机组，6条主要线路，表5-13和表5-14中的方案0给出了各机组的当前出力和各线路上对应的有功潮流值，方案1~32给出了围绕方案0的一些实验数据，试用这些数据确定各线路上有功潮流关于各发电机组出力的近似表达式。 2.设计一种简明、合理的阻塞费用计算规则，除考虑上述电力市场规则外，还需注意：在输电阻塞发生时公平地对待序内容量不能出力的部分和报价高于清算价的序外容量出力的部分。 3.假设下一个时段预报的负荷需求是982.4MW，表5-15、表5-16和表5-17分别给出了各机组的段容量、段价和爬坡速率的数据，试按照电力市场规则给出下一个时段各机组的出力分配预案。

第5章综合性实验 4.按照表5-18给出的潮流限值，检查得到的出力分配预案是否会引起输电阻塞，并在发生输电阻塞时，根据安全且经济的原则，调整各机组出力分配方案，并给出与该方案相应的阻塞费用。 5.假设下一个时段预报的负荷需求是1052.8MW，重复3~4的工作。

表5-13 各机组出力方案（单位：兆瓦，记作MW）
第5章综合性实验表5-13 各机组出力方案（单位：兆瓦，记作MW）方案\机组 1 2 3 4 5 6 7 8 120 73 180 80 125 81.1 90 133.02 129.63 158.77 145.32 78.596 75.45 … 28 75.529 29 104.84 30 111.22 31 98.092 32 120.44

表5-14 各线路的潮流值（各方案与表1相对应，单位：MW）
第5章综合性实验表5-14 各线路的潮流值（各方案与表1相对应，单位：MW）方案\线路 1 2 3 4 5 6 164.78 140.87 119.09 135.44 157.69 165.81 140.13 118.63 135.37 160.76 165.51 140.25 118.7 135.33 159.98 167.93 138.71 117.72 135.41 166.81 166.79 139.45 118.13 163.64 164.94 141.5 118.43 136.72 157.22 164.8 141.13 118.82 136.02 157.5 … 28 164.06 140.94 118.24 135.4 156.68 29 164.66 142.27 -147.2 120.21 135.28 157.65 30 164.7 142.94 120.68 135.16 157.63 31 164.67 141.56 119.68 135.29 157.61 32 164.69 143.84 121.34 135.12 157.64

第5章综合性实验表5-15 各机组的段容量（单位：MW）机组\段 1 2 3 4 5 6 7 8 9 10 70 50 30 40
50 30 40 20 15 110 55 75 95

表5-16 各机组的段价（单位：元/兆瓦小时，记作元/MWh）
第5章综合性实验表5-16 各机组的段价（单位：元/兆瓦小时，记作元/MWh）机组\段 1 2 3 4 5 6 7 8 9 10 -505 124 168 210 252 312 330 363 489 -560 182 203 245 300 320 360 410 495 -610 152 189 233 258 308 356 415 500 -500 150 170 200 255 302 325 380 435 800 -590 116 146 188 215 250 310 396 510 -607 159 173 205 305 405 520 120 180 251 260 306 315 335 348 548 -800 153 183 253 283 303 318 400

第5章综合性实验表5-17 各机组的爬坡速率（单位：MW/分钟）表5-18 各线路的潮流限值（单位：MW）和相对安全裕度机组 1
2 3 4 5 6 7 8 速率 2.2 3.2 1.3 1.8 1.4 线路 1 2 3 4 5 6 限值 165 150 160 155 132 162 安全裕度 13% 18% 9% 11% 15% 14%

第5章综合性实验【问题的假设】 1.在阻塞费用的计算中，仅考虑不能执行初始交易结果而付出的代价，而不考虑因输电阻塞引起的供电网络和其他硬件方面的损失； 2.各机组均不能低于其最低容量发电； 3.在初始交易结果确定后，无论阻塞管理怎样进行，原来合同内确定的各机组出力调整后仍然使用的容量部分均按该清算价格结算。

第5章综合性实验模型的建立线路潮流限制范围内的机组最大负荷和线路潮流裕度范围内的机组最大负荷的计算
(a)建立线路潮流限制范围内的机组最大负荷模型如下：目标函数：其中：式(3)为爬坡速率约束，式(4)为潮流限值约束

第5章综合性实验 (b)建立线路潮流裕度范围内的机组最大负荷模型如下：目标函数：
其中，式(5)为爬坡速率约束，式(6)为潮流限值裕度约束使用Matlab 软件求解以上两模型得：限值范围内极值为：；裕度范围内极值为：1094.6。

第5章综合性实验阻塞管理调整的模型建立 (a)第一步调整模型目标函数：阻塞费用其中：式(7)为负荷需求约束，
式(8)为爬坡速率约束，式(9)为潮流限值约束。

第5章综合性实验 (b)第二步调整模型在该步调整中，目标函数除阻塞费用外，还要考虑各线路上潮流值的绝对值超过限值的百分比。我们定义限值超过比来衡量。当时即第i 线路没有超过潮流限值时，限值超过比为0；当为时即第i 线路没有超过潮流限值时，限值超过比为，即：

第5章综合性实验目标函数：其中：式(10)为负荷需求约束，式(11)为爬坡速率约束，式(12)为线路潮流限值的裕度范围约束。

第5章综合性实验采用权和法将双目标非线性规划问题转化为单目标非线性规划问题如下：

第5章综合性实验第三步调整模型由于负荷需求超过了裕度范围内极值，为保证安全就必须拉闸限电。在裕度范围内最大负荷为1094.6，其相应的机组出力分配方案为：为保证最大限度的满足负荷需求，我们将裕度范围内的最大负荷作为负荷需求，将此方案作为初始交易结果，在负荷需求、爬坡速率、潮流裕度范围的约束下以阻塞费用为目标，可见此即第二步调整模型。由于此处已经拉闸限电，不能再使用简化的阻塞费用计算公式，具体模型如下：机组 1 2 3 4 5 6 7 8 出力 153 88 228 99.5 152 155 102.1 117

第5章综合性实验由问题三中的初始交易方案计算得各线路潮流值如下表其中，，为初始交易购电费用。线路 1 2 3 4 5 6 潮流值
其中，，为初始交易购电费用。由问题三中的初始交易方案计算得各线路潮流值如下表线路 1 2 3 4 5 6 潮流值

第5章综合性实验由上表与表6 比较可知有某些线路潮流超过了限值但没有超过裕度，因此必须进行阻塞管理。本题中负荷需求982.4 大于限值范围内极值，因此要进行第二步开始调整。我们通过调整λ值使用Matlab求解第二步调整模型结果如下：

第5章综合性实验当λ=200 时，阻塞费用为681.2，相应的调整后机组出力分配方案为：
当λ=20 时，阻塞费用为468.8，相应的调整后机组出力分配方案为：当λ=200 时，阻塞费用为681.2，相应的调整后机组出力分配方案为：当λ=2000 时，阻塞费用为1524.9，相应的调整后机组出力分配方案为：机组 1 2 3 4 5 6 7 出力（MW） 150 79 180 90 125 146 95 机组 1 2 3 4 5 6 7 出力（MW） 150 79 180 90 125 98.4 机组 1 2 3 4 5 6 7 出力（MW） 150 81 125 102.1

第5章综合性实验由计算结果，可见随权值λ增大，阻塞费用越大。这就表明对安全性重视程度变大后，阻塞费用就会相应的增加。

Similar presentations

Similar presentations

About project

反馈

请登录

Auth with social network:

Similar presentations

Similar presentations

About project

反馈