大连理工大学-矩阵大作业
- 格式:doc
- 大小:702.50 KB
- 文档页数:25
第三章 逐次逼近法1.1内容提要1、一元迭代法x n+1=φ(x n )收敛条件为:1)映内性x ∈[a,b],φ(x) ∈[a,b] 2)压缩性∣φ(x) -φ(y)∣≤L ∣x-y ∣其中L <1,此时φ为压缩算子,在不断的迭代中,就可以得到最终的不动点集。
由微分中值定理,如果∣φ’∣≤L <1,显然它一定满足压缩性条件。
2、多元迭代法x n+1=φ(x n )收敛条件为:1)映内性x n ∈Ω,φ(x n ) ∈Ω 2)压缩性ρ(▽φ)<1,其中▽φ为x n 处的梯度矩阵,此时φ为压缩算子,在不断的迭代中,就可以得到最终的不动点集。
3、当φ(x )= Bx+f 时,收敛条件为,ρ(B )<1,此时x n+1= Bx n +f ,在不断的迭代中,就可以得到线性方程组的解。
4、线性方程组的迭代解法,先作矩阵变换 U L D A --= Jacobi 迭代公式的矩阵形式 f Bx b D x U L D x n n n +=++=--+111)(Gauss-Seidel 迭代公式的矩阵形式 f Bx b L D Ux L D x n n n +=-+-=--+111)()( 超松弛迭代法公式的矩阵形式f Bx b L D x U D L D x k k k +=-++--=--+ωωωωω111)(])1[()(三种迭代方法当1)(<B ρ时都收敛。
5、线性方程组的迭代解法,如果A 严格对角占优,则Jacob 法和Gauss-Seidel 法都收敛。
6、线性方程组的迭代解法,如果A 不可约对角占优,则Gauss-Seidel 法收敛。
7、Newton 迭代法,单根为二阶收敛 2211'''21lim)(2)(lim---∞→+∞→--=-==--k k k k k k k k x x x x f f c x x ξξαα8、Newton 法迭代时,遇到重根,迭代变成线性收敛,如果知道重数m , )()('1k k k k x f x f m x x -=+仍为二阶收敛 9、弦割法)()())((111--+---=k k k k k k k x f x f x x x f x x 的收敛阶为1.618,分半法的收敛速度为(b-a )/2n-110、Aitken 加速公式11211112)(),(),(+----+-+--+---+---===k k k k k k k k k k k x x x x x x x x x x x ϕϕ1.2 典型例题分析1、证明如果A 严格对角占优,则Jacob 法和Gauss-Seidel 法都收敛。
作业1. 线性表编程作业:1.将顺序表逆置,要求用最少的附加空间。
参考答案#include <>#include <>#include <>#define LIST_INIT_SIZE 100#define LISTINCREMENT 10#define TRUE 1#define FALSE 0#define OK 1#define ERROR 0#define INFEASIBLE -1#define OVERFLOW -2typedef int Status;typedef int ElemType;typedef struct{ ElemType *elem;int length;int listsize;}SqList;立单链表 ");printf("2.取元素值 ");printf("3.查找 \n");printf("4.插入 ");printf("5.删除 ");printf("6.显示\n");printf("7.删除大于mink且小于maxk的元素值 ");printf("8.就地升序排序\n");printf("9.就地逆置 ");printf("a.有序表插入 ");printf("q.退出\n");printf("\n请选择操作:");fflush(stdin);scanf("%c",&choice);switch(choice){case '1': printf("请输入单链表中结点个数:");scanf("%d",&n);Create_L2(L,n);break;case '2': printf("请输入元素位序:");scanf("%d",&i);GetElem_L(L,i,e);printf("元素值为:%d\n",e);break;case '3': printf("请输入要查找的元素:");scanf("%d",&e);if(dlbcz(L,e))printf("查找成功!");elseprintf("查找失败。
第二章 矩阵变换和计算一、内容提要本章以矩阵的各种分解变换为主要内容,介绍数值线性代数中的两个基本问题:线性方程组的求解和特征系统的计算,属于算法中的直接法。
基本思想为将计算复杂的一般矩阵分解为较容易计算的三角形矩阵. 要求掌握Gauss (列主元)消去法、矩阵的(带列主元的)LU 分解、平方根法、追赶法、条件数与误差分析、QR 分解、Shur 分解、Jordan 分解和奇异值分解.(一) 矩阵的三角分解及其应用 1.矩阵的三角分解及其应用考虑一个n 阶线性方程组b Ax =的求解,当系数矩阵具有如下三种特殊形状:对角矩阵D ,下三角矩阵L 和上三角矩阵U ,这时方程的求解将会变得简单. ⎪⎪⎪⎪⎪⎭⎫⎝⎛=n d dd D21, ⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=nnn n l l l l l l L21222111, ⎪⎪⎪⎪⎪⎭⎫⎝⎛=nn n n u u u u u u U22212111. 对于b Dx =,可得解为i i i d b x /=,n i ,,2,1 =.对于b Lx =,可得解为1111/l b x =,ii i k k iki i l x lb x /)(11∑-=-=,n i ,,3,2 =.对于b Ux =,可得解为nn n n l b x /=,ii ni k k iki i l x lb x /)(1∑+=-=,1,,2,1 --=n n i .虽然对角矩阵的计算最为简单,但是过于特殊,任意非奇异矩阵并不都能对角化,因此较为普适的方法是对矩阵进行三角分解.1).Gauss 消去法只通过一系列的初等行变换将增广矩阵)|(b A 化成上三角矩阵)|(c U ,然后通过回代求与b Ax =同解的上三角方程组c Ux =的解.其中第k 步消元过程中,在第1-k 步得到的矩阵)1(-k A 的主对角元素)1(-k kka 称为主元.从)1(-k A 的第j 行减去第k 行的倍数)1()1(--=k kkk jkjk a a l (n j k ≤<)称为行乘数(子).2).矩阵A 的LU 分解对于n 阶方阵A ,如果存在n 阶单位下三角矩阵L 和n 阶上三角矩阵U ,使得LU A =, 则称其为矩阵A 的LU 分解,也称为Doolittle 分解.Gauss 消去法对应的矩阵形式即为LU 分解, 其中L 为所有行乘子组成的单位下三角矩阵, U 为Gauss 消去法结束后得到的上三角矩阵. 原方程组b Ax =分解为两个三角形方程组⎩⎨⎧==yUx b Ly .3).矩阵LU 分解的的存在和唯一性如果n 阶矩阵A 的各阶顺序主子式),,2,1(n k k =D 均不为零, 则必有单位下三角矩阵L 和上三角矩阵U ,使得LU A =, 而且L 和U 是唯一存在的.4).Gauss 列主元消去法矩阵每一列主对角元以下(含主对角元)的元素中, 绝对值最大的数称为列主元. 为避免小主元作除数、或0作分母,在消元过程中,每一步都按列选主元的Guass 消去法称为Gauss 列主元消去法.由于选取列主元使得每一个行乘子均为模不超过1的数,因此它避免了出现大的行乘子而引起的有效数字的损失.5).带列主元的LU 分解Gauss 列主元消去法对应的矩阵形式即为带列主元的LU 分解,选主元的过程即为矩阵的行置换. 因此, 对任意n 阶矩阵A ,均存在置换矩阵P 、单位下三角矩阵L 和上三角矩阵U ,使得LU PA =.由于选列主元的方式不唯一, 因此置换矩阵P 也是不唯一的. 原方程组b Ax =两边同时乘以矩阵P 得到Pb PAx =, 再分解为两个三角形方程组⎩⎨⎧==y Ux PbLy .5).平方根法(对称矩阵的Cholesky 分解)对任意n 阶对称正定矩阵A ,均存在下三角矩阵L 使T LL A =,称其为对称正定矩阵A 的Cholesky 分解. 进一步地, 如果规定L 的对角元为正数,则L 是唯一确定的.原方程组b Ax =分解为两个三角形方程组⎩⎨⎧==y x L bLy T .利用矩阵乘法规则和L 的下三角结构可得21112⎪⎪⎭⎫ ⎝⎛-=∑-=j k jkjj jjla l , jj j k jk ik ij ij l l l a l /11⎪⎪⎭⎫⎝⎛-=∑-=, i=j +1, j +2,…,n , j =1,2,…,n . 计算次序为nn n n l l l l l l l ,,,,,,,,,2322212111 .由于jj jk a l ≤,k =1,2,…,j .因此在分解过程中L 的元素的数量级不会增长,故平方根法通常是数值稳定的,不必选主元.6).求解三对角矩阵的追赶法 对于三对角矩阵⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛=---n nn n n b a c b a c b a c b 11122211A , 它的LU 分解可以得到两个只有两条对角元素非零的三角形矩阵 ⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛=⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛=--n n n nu d u d u d u l l l 11221132,1111U L . 其中⎪⎪⎩⎪⎪⎨⎧=-====-==--n i c l b u n i u a l b u n i c d i i i i i i i i i ,,3,2,,,3,2,/1,,2,1,1111计算次序是n n u l u l u l u →→→→→→→ 33221. 原方程组b Ax =分解为两个三角形方程组⎩⎨⎧==y Ux b Ly . 计算公式为n i y l b y b y i i i i ,,3,2,,111 =-==-,.1,,2,1,/)(,/1 --=-==+n n i u x c y x u y x i i i i i nn n该计算公式称为求解三对角形方程组的追赶法.当A 严格对角占优时,方程组b Ax =可用追赶法求解, 解存在唯一且数值稳定.7).矩阵的条件数设A 为非奇异矩阵,⋅为矩阵的算子范数,称1)(cond -=A A A 为矩阵A 的条件数.矩阵的条件数是线性方程组b Ax =, 当A 或b 的元素发生微小变化,引起方程组解的变化的定量描述, 因此是刻画矩阵和方程组性态的量. 条件数越大, 矩阵和方程组越为病态, 反之越小为良态.常用的矩阵条件数为∞-条件数: ∞-∞∞=1)(cond AA A ,1-条件数: 1111)(cond -=AAA ,2-条件数: )()()(cond mi n max 2122A A A A AAA HHλλ==-.矩阵的条件数具有如下的性质: (1) 1)(cond ≥A ;(2) )(cond )(cond 1-=A A ;(3) )(cond )(cond A A =α,0≠α,R ∈α;(4) 如果U 为正交矩阵,则1)(cond 2=U ,)(cond )(cond )(cond 222A AU UA ==.一般情况下,系数矩阵和右端项的扰动对解的影响为定理 2.5 设b Ax =,A 为非奇异矩阵,b 为非零向量且A 和b 均有扰动.若A 的扰动δA 非常小,使得11<-A A δ,则)()(cond 1)(cond bδb AδA AA A A xδx +-≤δ.关于近似解的余量与它的相对误差间的关系有定理2.6 设b Ax =,A 为非奇异矩阵,b 为非零向量,则方程组近似解x ~的事后估计式为bx A b A xx x bx A b A ~)cond(~~)cond(1-≤-≤-.其中称x A b ~-为近似解x ~的余量,简称余量。
矩阵与数值分析学生:学号:任课老师:金光日教学班号:(2)班院系:电子信息与电气工程学部《矩阵与数值分析》课程数值实验题目1.给定n 阶方程组A x b =,其中6186186186A ⎛⎫ ⎪ ⎪⎪= ⎪ ⎪ ⎪⎝⎭,7151514b ⎛⎫ ⎪⎪ ⎪= ⎪ ⎪⎪⎝⎭则方程组有解(1,1,,1)T x = 。
对10n =和84n =,分别用Gauss 消去法和列主元消去法解方程组,并比较计算结果。
1答: 程序1. Gauss 消元法function x=DelGauss(A,b) % Gauss 消去法 [n,m]=size(A); det=1; %存储行列式值 x=zeros(n,1); for k=1:n-1 for i=k+1:n if A(k,k)==0 return endm=A(i,k)/A(k,k); for j=k+1:nA(i,j)=A(i,j)-m*A(k,j); endb(i)=b(i)-m*b(k); enddet=det*A(k,k); %计算行列式enddet=det*A(n,n);for k=n:-1:1 %回代求解for j=k+1:nb(k)=b(k)-A(k,j)*x(j);endx(k)=b(k)/A(k,k);end2. 列主元Gauss消去法:function x=detGauss(A,b)% Gauss列主元消去法[n,m]=size(A);nb=length(b);det=1; %存储行列式值x=zeros(n,1);for k=1:n-1amax=0; %选主元for i=k:nif abs(A(i,k))>amaxamax=abs(A(i,k));r=i;endendif amax<1e-10return;endif r>k %交换两行for j=k:nz=A(k,j);A(k,j)=A(r,j);A(r,j)=z;endz=b(k);b(k)=b(r);b(r)=z;det=-det;endfor i=k+1:n %进行消元m=A(i,k)/A(k,k);for j=k+1:nA(i,j)=A(i,j)-m*A(k,j);endb(i)=b(i)-m*b(k);enddet=det*A(k,k);enddet=det*A(n,n);for k=n:-1:1 %回代求解for j=k+1:nb(k)=b(k)-A(k,j)*x(j);endx(k)=b(k)/A(k,k);end矩阵A和b的构造clc;clear;n=10;%n=84;A=eye(n)*6+diag(ones(1,n-1)*8,-1)+diag(ones(1,n-1),1); b=[7,15*ones(1,n-2),14]';计算结果:(1)n=10时Gauss消元法>>x=DelGauss(A,b)x =1.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000列主元Gauss消去法>>x=detGauss(A,b)x =1111111111(2) n=84时Gauss消元法>>x=DelGauss(A,b) x =1.0e+008 *0.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0001 0.0002 -0.0003 0.0007 -0.0013 0.0026 -0.0052 0.0105 -0.0209 0.0419 -0.0836 0.16650.6501-1.25822.3487-4.02635.3684列主元Gauss消去法>>x=detGauss(A,b) x =1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.00001.0000 1.0000 1.00001.00001.0000 1.0000结果分析由上述实验结果可知,对于n=10采用Gauss 消去法和Gauss 列主元消去法得到的实验结果是相同的,而对于n=84,Gauss 消去法所得到的结果是错误的,Gauss 列主元消去法得到的结果是正确的。
习题2-11. =6.32A 2. 用行列式的定义计算下面的行列式.(1)35;(2)256;(3)8;(4)29.−−思考题 2-21.若对方阵A 进行一次对调变换得到,则B =−A B ;若对方阵A 进行一次倍乘变换(假设第i 行或第i 列乘以数)得到,则k B k =B A ;若对方阵A 进行一次倍加变换得到,则B .=A B2.0.=A3.(1)不正确。
例如,设则 1112111221222122,,a a b b a a b b ⎡⎤⎡==⎢⎥⎢⎣⎦⎣A B ⎤⎥⎦1111121211121211121221212222212222212222a b a b a a b b a b a b a b a a b b a b +++++==+++++A B111211121112111211121112212221222122212221222122a a ab b a b b a b b aa a ab b a b b a b b a =+++=+++A B(2)不正确。
设A 的阶数为,则n (1)n−=−A A (3)不正确。
例如,设,则1200⎡⎤=⎢⎣⎦A ⎥0,=A 但.≠A O 4. ,,1,(),()1i j i i j k k k =−==E E E5. 性质2-2讲的是方阵A 的第行(列)的数与第i 行(列)对应的代数余子式的乘积之和等于i A 的行列式;性质2-7讲的是方阵A 的第i 行(列)的数与另一行(列)对应的代数余子式的乘积之和等于0.习题2-21. 2111231123123det()3,,39,,9,,18.c c a a a a a a a a a a a −=+−=−+=−=−A 2. 131223123233122312312323,2,3,,3,,3,,c c c c c c −+−−++=−===a a a a a a a a a a a a a a a a 63.321123211321212311223,,,,,,,,,,,,,,,n m +=+=−+=−a a a b b a a a b a a a b a a a b a a b a4.证:(1)将第2列和第3列都加到第1列,得0000a b b c c a b c c ab c c a a b c a a b c a a b b ca b b c−−−−−−−−=−−=−−−−−. (2)111111111111111122222222222222223333333333333333a b b c c a a b c c a b b c c a a b b c c a a b c c a b b c c a a b b c c a a b c c a b b c c a ++++++++++=++++++++++++ 1111111111111111122222222222222222333333333333333332a b c c b c c a a b c b c a a b c a b c c b c c a a b c b c a a b c a b c c b c c a a b c b c a a b c ++=+++=+=++ (3)设A 的阶数为,则为奇数.由n n A 是反称矩阵,得T=−A A .两边取行列式,得 ,(1),Tn=−=−=−,A A A A A A 故0.=A 5. 先按行提公因式,在按列提公因式,得2111121211221212222221122n n n n n n n n nn na b a b b a b b a b b a b a b b a b b a b b a b11112212112222121122n n n nn n n nn a b a b a b a ba b a b b b b a b a b a b =n1112121222222222121212n nnn n n nna a a a a ab b bb b bc a a a ==6.(1)解:先按行提公因式,在按列提公因式,得1111114111ab ac ae bd cd de abcdef abcdef bfcfef −−−=−=−−(2)103100204310043141992003951200510012520301300600130013=−−=−−=提高题2-21.,,,,,,+=++++=+−++A B ξηαββγαγξηαγβγαγ,,,,,,22,,,=+−++=+−+=+ξηαγβγαγξηαγβγγξηαβγ2(,,,,,,)2()6=+=+ξαβγηαβγA B =2.1231231231232323,24,36,3,25=++++++=++++B a a a a a a a a a a a a a a a a 1232331223123,3,,,,,=+++−=−+=−=−a a a a a a a a a a a a a 103.根据性质2-7,得41424344414243441111A A A A A A A A +++=⋅+⋅+⋅+⋅=4.(1).132343(1)(1)52(1)301(1)415D +++=−⋅−+⋅−++⋅−=− (2) 1424449(1)(1)52(1)01(1)40,2a a +++−⋅−+⋅−++⋅−==−.5.(1)对第2行和第4行分别应用性质2-2和性质2-7,得212223242521222324254()3()4,2()()0A A A A A A A A A A ++++=⎧⎨++++=⎩ 解得.2122232A A A ++=−(2)对第2行和第4行分别应用性质2-7,得313233343531323334354()3()0,2()()0A A A A A A A A A A ++++=⎧⎨++++=⎩解得=0.313233A A A ++思考题 2-31.表示第二行先乘以2,再用第二行减去第一行,22r r −12122323112012r r −=.2.对行列式进行对调变换和倍乘变换时,需要在得出的行列式的前面添加负号和系数,对行列式进行初等变换时,关心的是最后的数值;对矩阵进行初等变换时不需要添加负号和系数,对矩阵进行初等变换时,关心的是用何种变换进行化简,最后化成何种形式。
(1)从大到小的顺序的计算程序:function y=snd(n)format longy=0;if n<2disp('请输入大于1的数!')elses=0;i=2;while i<=ns=single(s+(1/(i^2-1))); i=i+1;endy=s;end(2)从小到大的顺序的计算程序:function y=snx(n)format longy=0;if n<2disp('请输入大于1的数!')elses=0;i=n;while 1s=single(s+(1/(i^2-1)));i=i-1;if i==1breakendendy=s;end(3)按两种顺序分别计算并指出有效位数(编制程序时用单精度)S的计算结果:①210S的计算结果:②410S的计算结果:③610计算时的有效位数为七位数。
① 秦九昭算法计算程序:function y=qjz(a,x) j=3; i=size(a,2); switch i case 1 y=a(1); case 2y=a(1)*x+a(2); otherwise p=a(1)*x+a(2); while j<=i p=p*x+a(j); j=j+1; end y=p; end② 计算在点23的值。
计算结果如下:当23=x 时()86652=x f 。
①Gauss法计算程序和结果:程序:A(1,:)=[31,-13,0,0,0,-10,0,0,0]; A(2,:)=[-13,35,-9,0,-11,0,0,0,0]; A(3,:)=[0,-9,31,-10,0,0,0,0,0];A(4,:)=[0,0,-10,79,-30,0,0,0,-9]; A(5,:)=[0,0,0,-30,57,-7,0,-5,0]; A(6,:)=[0,0,0,0,-7,47,-30,0,0];A(7,:)=[0,0,0,0,0,-30,41,0,0];A(8,:)=[0,0,0,0,-5,0,0,27,-2];A(9,:)=[0,0,0,-9,0,0,0,-2,29];B=[-15;27;-23;0;-20;12;-7;7;10]; [a,b]=gauss(A,B);j=size(a,2);while j>=1k=j+1;s=b(j);while k<=9s=s-x(k)*a(j,k);k=k+1;endx(j)=s/a(j,j);j=j-1;enddisp(x)function [x,y]=gauss(a,b)num_i=size(a,1);j=1;while j<=(num_i-1)i=j+1;while i<=num_ir=a(i,j)/a(j,j);a(i,:)=a(i,:)-r*a(j,:);b(i,:)=b(i,:)-r*b(j,:);i=i+1;endj=j+1;endx=a;y=b;运行的结果为:()T 289.0-345.0.0=。
(1)从大到小的顺序的计算程序:function y=snd(n)format longy=0;if n<2disp('请输入大于1的数!')elses=0;i=2;while i<=ns=single(s+(1/(i^2-1))); i=i+1;endy=s;end(2)从小到大的顺序的计算程序:function y=snx(n)format longy=0;if n<2disp('请输入大于1的数!')elses=0;i=n;while 1s=single(s+(1/(i^2-1)));i=i-1;if i==1breakendendy=s;end(3)按两种顺序分别计算并指出有效位数(编制程序时用单精度)S的计算结果:①210S的计算结果:②410S的计算结果:③610计算时的有效位数为七位数。
① 秦九昭算法计算程序:function y=qjz(a,x) j=3;i=size(a,2); switch i case 1 y=a(1); case 2y=a(1)*x+a(2); otherwisep=a(1)*x+a(2); while j<=i p=p*x+a(j); j=j+1; end y=p; end② 计算在点23的值。
计算结果如下:当23=x 时()86652=x f 。
①Gauss法计算程序和结果:程序:A(1,:)=[31,-13,0,0,0,-10,0,0,0]; A(2,:)=[-13,35,-9,0,-11,0,0,0,0]; A(3,:)=[0,-9,31,-10,0,0,0,0,0];A(4,:)=[0,0,-10,79,-30,0,0,0,-9]; A(5,:)=[0,0,0,-30,57,-7,0,-5,0]; A(6,:)=[0,0,0,0,-7,47,-30,0,0];A(7,:)=[0,0,0,0,0,-30,41,0,0];A(8,:)=[0,0,0,0,-5,0,0,27,-2];A(9,:)=[0,0,0,-9,0,0,0,-2,29];B=[-15;27;-23;0;-20;12;-7;7;10]; [a,b]=gauss(A,B);j=size(a,2);while j>=1k=j+1;s=b(j);while k<=9s=s-x(k)*a(j,k);k=k+1;endx(j)=s/a(j,j);j=j-1;enddisp(x)function [x,y]=gauss(a,b)num_i=size(a,1);j=1;while j<=(num_i-1)i=j+1;while i<=num_ir=a(i,j)/a(j,j);a(i,:)=a(i,:)-r*a(j,:);b(i,:)=b(i,:)-r*b(j,:);i=i+1;endj=j+1;endx=a;y=b;运行的结果为:()T 289.0-345.0.0=。
713----154.0201.0x290058.0221.0430.0.0②列主元消去法计算程序和结果:A(1,:)=[31,-13,0,0,0,-10,0,0,0];A(2,:)=[-13,35,-9,0,-11,0,0,0,0];A(3,:)=[0,-9,31,-10,0,0,0,0,0];A(4,:)=[0,0,-10,79,-30,0,0,0,-9];A(5,:)=[0,0,0,-30,57,-7,0,-5,0];A(6,:)=[0,0,0,0,-7,47,-30,0,0];A(7,:)=[0,0,0,0,0,-30,41,0,0];A(8,:)=[0,0,0,0,-5,0,0,27,-2];A(9,:)=[0,0,0,-9,0,0,0,-2,29];B=[-15;27;-23;0;-20;12;-7;7;10];[a,b]=lzy(A,B);j=size(a,2);while j>=1k=j+1;s=b(j);while k<=9s=s-x(k)*a(j,k);k=k+1;endx(j)=s/a(j,j);j=j-1;enddisp(x)function [a1,b1]=lzy(a,b)[num_i,num_j]=size(a);ab=zeros(num_i,num_j+1);for k=1:num_jab(:,k)=a(:,k);endab(:,num_j+1)=b(:,1);j=1;while j<num_j[max,max_i]=searmax(j,j,ab);i_nu=ab(j,:);ab(j,:)=ab(max_i,:);ab(max_i,:)=i_nu;m=j+1;while m<=num_ifor n=j:num_j+1ab(m,n)=ab(m,n)-(ab(m,j)/max)*ab(j,n);endm=m+1;endj=j+1;enda1=zeros(num_i,num_j);for k=1:num_ia1(:,k)=ab(:,k);endb1=ab(:,num_i+1);function [b,c]=searmax(i,j,a)num_i=size(a,1);k=i;m=abs(a(k,j));c=k;while k<num_ik=k+1;if m>=abs(a(k,j))continueelsem=abs(a(k,j));c=k;endendb=a(c,j);运行的结果为:()T 477.0767237.0-=---.0-.0078285.0x345171.0.0308.0146.0(1)LU分解的计算程序及结果:function [l,u]=lufz(a)[num_i,num_j]=size(a);l=eye(num_i); u=eye(num_i); for j=1:num_j u(1,j)=a(1,j); l(j,1)=a(j,1)/u(1,1); end i=2;while i<=num_i j=i;while j<num_i s=0; for k=1:i-1s=s+l(i,k)*u(k,j); endu(i,j)=a(i,j)-s; s=0; for k=1:i-1s=s+l(j+1,k)*u(k,i); endl(j+1,i)=(a(j+1,i)-s)/u(i,i); j=j+1; end s=0; for k=1:i-1s=s+l(i,k)*u(k,num_i); endu(i,num_i)=a(i,num_i)-s; i=i+1; end输入要求解的矩阵后运行的结果为:⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛--------------=1092.0020.0014.0083.0119.00001025.0018.0112.000001655.0000001157.000001398.00001354.0001305.001419.01L⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭⎝------------------=390.27420.2413.26367.0513.0381.21562.0785.0308732.45578.350180.7602.44900452.0186.31461.75000277.1350.310259.28000194.41109548.29U(2)带列主元的LU 分解计算程序及结果由于Matlab 中自带带列主元的LU 分解函数,故计算程序如下:A(1,:)=[31,-13,0,0,0,-10,0,0,0]; A(2,:)=[-13,35,-9,0,-11,0,0,0,0]; A(3,:)=[0,-9,31,-10,0,0,0,0,0]; A(4,:)=[0,0,-10,79,-30,0,0,0,-9]; A(5,:)=[0,0,0,-30,57,-7,0,-5,0]; A(6,:)=[0,0,0,0,-7,47,-30,0,0]; A(7,:)=[0,0,0,0,0,-30,41,0,0]; A(8,:)=[0,0,0,0,-5,0,0,27,-2]; A(9,:)=[0,0,0,-9,0,0,0,-2,29]; [L,U,P]=lu(A);运行结果如下:⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛--------------=1092.0020.0014.0083.0119.00001025.0018.0112.000001655.0000001157.000001398.00001354.0001305.001419.01L⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭⎝------------------=390.27420.2413.26367.0513.0381.21562.0785.0308732.45578.350180.7602.44900452.0186.31461.75000277.1350.310259.28000194.41109548.29U ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛=100000000010000000001000000000100000000010000000001000000000100000000010000000001P 为单位阵。
(3)求A 的逆矩阵的计算程序及结果:A(1,:)=[31,-13,0,0,0,-10,0,0,0]; A(2,:)=[-13,35,-9,0,-11,0,0,0,0]; A(3,:)=[0,-9,31,-10,0,0,0,0,0]; A(4,:)=[0,0,-10,79,-30,0,0,0,-9]; A(5,:)=[0,0,0,-30,57,-7,0,-5,0]; A(6,:)=[0,0,0,0,-7,47,-30,0,0]; A(7,:)=[0,0,0,0,0,-30,41,0,0]; A(8,:)=[0,0,0,0,-5,0,0,27,-2]; A(9,:)=[0,0,0,-9,0,0,0,-2,29]; [num_i,num_j]=size(A); A_n=eye(num_i); E=eye(num_i); [L,U]=lufz(A); for i=1:num_i y=hd(1,L,E(:,i)); A_n(:,i)=hd(0,U,y); end disp(A_n)function [l,u]=lufz(a) [num_i,num_j]=size(a);l=eye(num_i);u=eye(num_i);for j=1:num_ju(1,j)=a(1,j);l(j,1)=a(j,1)/u(1,1);endi=2;while i<=num_ij=i;while j<num_is=0;for k=1:i-1s=s+l(i,k)*u(k,j);endu(i,j)=a(i,j)-s;s=0;for k=1:i-1s=s+l(j+1,k)*u(k,i);endl(j+1,i)=(a(j+1,i)-s)/u(i,i);j=j+1;ends=0;for k=1:i-1s=s+l(i,k)*u(k,num_i);endu(i,num_i)=a(i,num_i)-s;i=i+1;endfunction x=hd(f,a,b)[num_i,num_j]=size(a);x=zeros(num_i,1);switch fcase 0i=num_i-1;x(num_i)=b(num_i)/a(num_i,num_j); while i>=1s=0;for k=i+1:num_is=s+a(i,k)*x(k);endx(i)=(b(i)-s)/a(i,i);i=i-1;end case 1x(1)=b(1)/a(1,1); j=2;while j<=num_i s=0; for k=1:j-1s=s+a(j,k)*x(k); endx(j)=(b(j)-s)/a(j,j); j=j+1; end otherwisedisp('error !请输入正确的参数!') end运行结果:⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛=-0365.00034.00008.00011.00036.00058.00020.00006.00003.00033.00382.00010.00014.00048.00023.00008.00002.00001.00007.00010.00468.00306.00050.00021.00007.00002.000010.00013.00306.00419.00068.00028.00010.00003.00001.00035.00048.00051.00070.00243.00101.00036.00011.00005.00058.00024.00024.00033.00105.00181.00064.00020.00008.00025.00015.00028.00039.00069.00077.00381.00116.00049.00022.00024.00071.00097.00121.00065.00131.00378.00159.00012.00014.00129.00176.00073.00036.00058.00160.00390.01A(4)求A 的行列式的计算程序及结果:A(1,:)=[31,-13,0,0,0,-10,0,0,0]; A(2,:)=[-13,35,-9,0,-11,0,0,0,0]; A(3,:)=[0,-9,31,-10,0,0,0,0,0]; A(4,:)=[0,0,-10,79,-30,0,0,0,-9]; A(5,:)=[0,0,0,-30,57,-7,0,-5,0]; A(6,:)=[0,0,0,0,-7,47,-30,0,0]; A(7,:)=[0,0,0,0,0,-30,41,0,0]; A(8,:)=[0,0,0,0,-5,0,0,27,-2]; A(9,:)=[0,0,0,-9,0,0,0,-2,29]; [num_i,num_j]=size(A); [L,U]=lufz(A); s=1;for i=1:num_i s=s*U(i,i); enddisp(['矩阵的行列式值为:',num2str(s)])运行程序后结果与调用matlab中det()函数结果如下:(1)求cholesky分解程序及结果:function l=chole(a)[num_i,num_j]=size(a);l=zeros(num_i);l(1,1)=sqrt(a(1,1));for k=2:num_il(k,1)=a(k,1)/l(1,1);endj=2;while j<=num_js1=0;for k=1:j-1s1=s1+l(j,k)^2;endl(j,j)=sqrt(a(j,j)-s1);i=j+1;while i<=num_is2=0;for k=1:j-1s2=s2+l(i,k)*l(j,k);endl(i,j)=(a(i,j)-s2)/l(j,j);i=i+1;endj=j+1;end运行程序后的结果为:⎪⎪⎪⎪⎪⎭⎫ ⎝⎛--=7494.07413.00641.37071.008480.21785.17071.0001213.27071.00004142.1L (2)求解方程组程序及结果:a=[2,1,-1,1;1,5,2,7;-1,2,10,1;1,7,1,11]; b=[13;-9;6;0]; L=chole(a); y=hd(1,L,b); x=hd(0,L',y); disp(x)运行后的结果为:()0732.385122.123902.654146.26-=x程序和结果如下:A=[4,2,0,0;3,-2,1,0;0,2,5,3;0,0,-1,6]; B=[6;2;10;5]; [L,U]=sjfj(A);% y=hd(1,L,B); x=hd(0,U,y); disp('方程组的解为:')disp(x)function [l,u]=sjfj(a) num_i=size(a,1); i=2;while i<=num_ia(i,i-1)=a(i,i-1)/a(i-1,i-1); a(i,i)=a(i,i)-a(i,i-1)*a(i-1,i); i=i+1; end l=tril(a); u=triu(a); for j=1:num_i; l(j,j)=1; end运行程序后结果为:()Tx 1111=编程求解矩阵A 的QR 分解:(1)QR 分解计算程序:function [Q,R]=hqr(A) [n,n]=size(A); Q=eye(n); for i=1:(n-1) E=eye(n-i+1); e1=E(:,1);a=zeros(1,n-i+1)'; for j=1:(n-i+1) a(j)=A(j+i-1,i); endav=sqrt(a'*a); w=a-av*e1;h=E-(2/(w'*w))*(w*w'); q=eye(n); for k=i:n for j=i:nq(k,j)=h(k-i+1,j-i+1); end end A=q*A; Q=q*Q; end R=A; Q=Q';(2)输入矩阵A 后的计算结果:⎪⎪⎪⎪⎪⎭⎫ ⎝⎛-------=5774.05164.006325.007746.006325.05774.02582.07071.03162.05774.02582.07071.03162.0Q ⎪⎪⎪⎪⎪⎭⎫ ⎝⎛------=1547.10000328.15492.1003536.03536.08284.200555.24743.01623.31623.3R(1)Jacobi 迭代法求解程序与结果:a=0; b=0; c=0; while 1a0=0.25*(7+b-c); b0=(1/8)*(-21-4*a-c); c0=(1/5)*(15+2*a-b);A=[abs(a-a0),abs(b-b0),abs(c-c0)]; f=max(A); if f<=0.001 x=[a0,b0,c0]; break else a=a0; b=b0; c=c0; end enddisp('方程组的解为x') disp(x)运行后的结果为:()Tx 647.3111.3060.0-=(2)Gauss-Seidel 迭代法求解的程序与结果:a=0; b=0; c=0; while 1 a0=a; b0=b; c0=c;a=(1/4)*(7+b-c); b=(1/8)*(-21-4*a-c); c=(1/5)*(15+2*a-b);A=[abs(a-a0),abs(b-b0),abs(c-c0)]; f=max(A); if f<=0.001 x=[a,b,c]; breakend enddisp('方程组的解为:') disp(x)运行后的结果为:()Tx 646.3111.3061.0-=(1)计算01523=--x x 在[]2,1上的根的程序:syms x f=2*x^3-5*x-1; df=diff(f,x); g=x-(f/df); x0=1; while 1x1=subs(g,x0); dx=abs(x1-x0); if dx<=0.001disp(['Newton 法的解为:x=',num2str(x1)]) break else x0=x1; end end x0=1; x1=1.1; while 1x2=x1-(subs(f,x1)/(subs(f,x1)-subs(f,x0)))*(x1-x0); dx=abs(x2-x1); if dx<=0.001disp(['割线法的解为:x=',num2str(x2)]) break else x0=x1; x1=x2; end end运行后结果为:Newton 法673.1=x ,割线法673.1=x 。