计算方法实验

1、数值计算方法实验专业:姓名:姓名:实验一1. 题目 (课本 46页习题二第 1 题) 用列主元素法解下列方程，结果保留 4 位小数2x1 6x2 4x34x1 4x2 5x3 3 6x1 x2 18x322. 模型建立(理论分析)用方程组的增广矩阵a11a12a1nb1A,ba21a22a2nb2an1an2annbn表示，并在增广矩阵上进行计算。首先在矩阵的第一列中选取绝对值最大的元，通过行互换，将最大的元所在行换在第一行的位置，得到增广矩阵A1 ,b1 ，然后进行消元，得到矩阵11 a11 a1211 a1nb122 A ,b2 a2222 a2nb2an222 annbn然后在每次从消元

2、之前，都通过行变换将绝对值最大的元所在的行换到相应 行，再进行消元。如此经过 n-1 步，增广矩阵被化成上三角矩阵，最后回代过程求 解。这种算法就是列主元素法求解方程组。3. 算法(1) 输入 A=2,6,-4;1,4,-5;6,-1,18;b=4;3;2; 计算 rank(A);rank(B) 。(2) 判断所需求的方程组解的情况，选择相应计算方法，本题是存在唯一解。(3) 置储存空间X, C矩阵，比较所在列进行消元运算数值绝对值，取下标位置。(4) 取绝对值大的元所在行通过行变换换到首位，再进行消元，直到化为上三角矩 阵，停机。4. 程序列主元素消去法程序 clearclcA=2,6,-4

3、;1,4,-5;6,-1,18; b=4;3;2;B=A b;RA=rank(A);RB=rank(B); zhica=RB-RA;if zhica0,disp(请注意：因为RA=RB，所以此方程组无解)returnendif RA=RBif RA=3disp(请注意：因为 RA=RB=3，所以此方程组有唯一解) X=zeros(3,1);C=zeros(1,4);for p=1:2 Y,j=max(abs(B(p:3,p); C=B(p,:);B(p,:)=B(j+p-1,:);B(j+p-1,:)=C;for k=p+1:3m=B(k,p)/B(p,p);B(k,p:4)=B(k,p:4)

4、-m*B(p,p:4);end end b=B(1:3,4);A=B(1:3,1:3);X(3)=b(3)/A(3,3);for q=2:-1:1X(q)=(b(q)-sum(A(q,q+1:3)*X(q+1:3)/A(q,q);endXelsedisp(请注意：因为 RA=RB3，所以此方程组有无穷多解)end end5. 运行结果祇 实验二1. 题目 （课本 46页习题二第 3 题）用紧凑格式解下列方程组，并写出 L，U 矩阵1234x1214916x210182764x34411681256x41902. 模型建立（理论分析）已知用矩阵形式表示 Gauss消去法的消元过程是对方程组增广矩

5、阵A,b进行一系列初等变换，将系数矩阵 A 左乘上三角形矩阵的过程，也等价于用一串初等 矩阵去做乘增广矩阵，因此消元过程可以通过矩阵运算来表示。例如，第一次消 元等价于用初等矩阵1000l21 100L1l31 01011左乘矩阵 A1，b 1A，bln1 001其中 li1ai11 /1 a11 i2,3,n ，即 A 2，b 2L1 A 1，b1 。般地，第 k 次消元等价于用初等矩阵101k 1k1l nk 0,n ，经过 n 1次消元后得到左乘矩阵 Ak，bk ，其中 lik aikk /akkk i k 1,lnk 011l211011LL11L21Ln11l31l321， A， b

6、nnLA n ，Lb n1ln1ln2l n3lnn 1 1An，bnLn 1Ln 2 L1 A1，b1 ， 因 为 Lk1 01lk 1k 1，令消元过程实际上是把系数矩阵 A 分解成单位下三角与上三角矩阵的乘积的过程。A LU其中 L 为单位上三角矩阵， U 为下三角矩阵。1u11u12u13u1nl21 10u22u23u2nl31 l32 1，UL 31 32ln1ln2ln31l nn 110unn这种分解称为杜利特尔(Doolittle) 分解，也称为LU分解。当A为n阶方阵，若A的顺序主子式A i 1,2, ,n 1均不为零，则矩阵A存在唯一的 Doolittle 分解。由矩阵的

7、乘法法则，得 lij和 uij的公式u1ja1j j 1,2,nuijaiji1likukj ik12,n, ji, ,nlijijaijj1likukj/ujj j1,2,n 1,i j 1, ,nk1计算过程应按第1行，第1列，第2行，第2列，的顺序进行。得到L , U矩阵 再由Ly b和Ux y计算方程组的解。本题要求用紧凑格式解下列方程组， 并写出 L， U 矩阵。在程序实现上就是进 行 Doolittle 分解。3. 算法(1) 输入 A=1,2,3,4;1,4,9,16;1,8,27,64;1,16,81,256;b=2;10;44;190;并设 l,u 的储存空 间。u1ja1j

8、 j 1,2, ,ni1(2) 计算 uijaijlikukj i 2,n, j i, ,nk1j1lijaijlikukj /ujj j 1,2, ,n 1,i j 1, ,nk1Ly b和Ux y计算方程组的解，停机4. 程序LU 矩阵分解clearclc A=1,2,3,4;1,4,9,16;1,8,27,64;1,16,81,256; b=2;10;44;190;n=length(A); u=zeros(n,n);l=eye(n,n); u(1, : )=A(1, : );l(2:n,1)=A(2:n,1)/u(1,1); for k=2:nu(k,k:n)=A(k,k:n)-l(k,

9、1:k-1)*u(1:k-1,k:n); l(k+1:n,k)=(A(k+1:n,k)-l(k+1:n,1:k-1)*u(1:k-1,k)/u(k,k); end ulY=zeros(n,1);Y(1)=b(1);for i=2:nY(i)=b(i)-l(i,1:i-1)*Y(1:i-1);end X=zeros(n,1); if det(u)=0;X=0;elseX(n)=Y(n)/u(n,n);for i=n-1:-1:1X(i)=(Y(i)-u(i,i+1:n)*X(i+1:n)/u(i,i);endend X5. 运行结果4- MATLAB 7.6.0 (RSsf *1=1回j-3l)

10、1Ale Edit Debu Parallel Desktop Window Help总町 4 誓& W -2014200093 - |_| 实验三1. 题目 （课本 47页习题二第 6 题）分别用平方根法和改进平方根法求解方程组1213x112505x2210141x31635115x482. 模型建立（理论分析）依据定理，如果A为n阶对称正定矩阵，则存在惟一的非奇异下三角阵L使得A LLT且L的对角元素皆为正数。1111 21n1A l21 l 22l 22ln2ln1 l n2l nnlnn其中 lii 0 i 1,2 n。由矩阵乘法及 l ik0（当 i k 时），得naijlik l

11、 jkk1j1k1lik l jk l jj lij矩阵的这种分法称为Cholesky 分解法，即平方根分解法。比较A与LLt的相应元素，可得k11/2l kkakklkjj1k1,2, ,nk1likaiklijlkjj1/lkki k 1, ,n规0定0。j1计算顺序是按列进行 的 ，即l11li1 i2,3, ,nl22li2 i 3, ,n当矩阵 A 完成 Cholesky 分解后，求解对称正定方程组Axb的解计算公式。ykbklkjyj,k 1,2, ,nj 1nXk(yk IjkXj) lkk,k n, n 1,1j k 1当A为对称矩阵，且A的所有顺序主子式均不为零， A可唯一分

12、解为如下的形式A LU为了利用A的对称性，将U再分解为:U11U22U12U11UnnU1nU11Un-1,nU e1,n-1DU。其中D为对角阵，Uo为单位上三角阵，于是 A LU LDU 0由分解的唯一性即得Uo LT，代入得到对称矩阵A的分解式A LDLT。这是一种 改进的平方根法（LDLT法）其中L为单位下三角阵，D diag d1, ,dn为对角阵。 由比较法可以导出LDLt分解的计算公式：对于k 1,2, ,nk 1dkakklj 12 d. kj jk 1hkaikj 1lij d jlkj/dk ik 1, ,n计算顺序如下：d1 li1 i 2,3,nd2li2i 3, ,n

13、da此时式中有许多计算式重复的,改进辅助量uik1hkdk k1,2,n,i k 1, ,nk 1dkakkUkjl j 1k 1kjUkjaikUijlkjj 1l ik ukj /dkLDLt分解后，解对称正定方程组 Ax b可分两步进行：先解方程组Ly b，再由LTx D 1 y 求 x。本问题使用平方根法与改进的平方根法两种方法解方程组。3. 算法两种解法算法应相似，只介绍平方根法k 1 1/2akkl k2jk1,2,j1k1aiklijlkj /lkkikj1lkk(2)计算lik(1)输入 A=1,2,1,-3;2,5,0,-5;1,0,14,1;-3,-5,1,15;b=1;2

14、;16;8;,n，得到 L 矩阵1, ,n(3) 求解先解方程组Ly b，再由Ltx4. 程序(1)平方根法求解方程组程序段clearclcA=1,2,1,-3;2,5,0,-5;1,0,14,1;-3,-5,1,15; b=1;2;16;8;n=length(b);L=zeros(n,n);L(1,1)=sqrt(A(1,1);L(2:n,1)=A(2:n,1)/L(1,1);for j=2:n-1L(j,j)=sqrt(A(j,j)-sum(L(j,1:j-1)A2);for i=j+1:nL(i,j)=(A(i,j)-sum(L(i,1:j-1).*L(j,1:j-1)/L(j,j);e

15、ndendL(n,n )=sqrt(A( n,n )-sum(L( n,1: n-1)42);Lx=zeros(n,1);y=zeros(n,1);y(1)=b(1)/L(1,1);for i=2:ny(i)=(b(i)-sum(L(i,1:i-1).*y(1:i-1)/L(i,i);endyL=L;x(n)=y(n)/L(n,n);for i=n-1:-1:1x(i)=(y(i)-sum(L(i,i+1:n).*x(i+1)/L(i,i);endX( 2)改进平方根法求解方程组程序段 clearclcA=1,2,1,-3;2,5,0,-5;1,0,14,1;-3,-5,1,15; b=1;2

16、;16;8;L=zeros(4,4);U=zeros(4,4);D=zeros(4,4);for i=1:4L(i,i)= 1;endD(1,1)=A(1,1);for i=2:4U(i,1)=A(i,1);L(i,1)= U(i,1)/D(1,1);endD(2,2)=A(2,2)-U(2,1)*L(2,1);for i=3:4U(i,2)=A(i,2)-U(i,1)*L(2,1);L(i,2)= U(i,2)/D(2,2);endD(3,3)=A(3,3)-U(3,1)*L(3,1)-U(3,2)*L(3,2);U(4,3)=A(4,3)-U(4,1)*L(3,1)-U(4,2)*L(3,

17、2);L(4,3)= U(4,3)/D(3,3);D(4,4)=A(4,4)-U(4,1)*L(4,1)-U(4,2)*L(4,2)-U(4,3)*L(4,3); LDx=zeros(4,1);y=zeros(4,1);y(1)=b(1)/L(1,1);for i=2:4 y(i)=(b(i)-sum(L(i,1:i-1).*y(1:i-1)/L(i,i);endyL=L;y=inv(D)*y;x(4)=y(4)/L(4,4);for i=3:-1:1x(i)=(y(i)-sum(L(i,i+1:4).*x(i+1)/L(i,i);endx5. 运行结果(1)平方根法求解方程组程序运行结果尸

18、亘 I实验四1.题目(课本72页习题三第4题)分别用Gauss-Seide迭代法与SOR法(s =1.25求解方程组4 30 xi2434-1 X2300 -14X324取x01 1 1T,迭代7次，并比较它们的计算结果2. 模型建立(理论分析)迭代法主要解决大型稀疏矩阵方程组，该方法包括Jacobi迭代法，Gauss-Seidel 迭代法，超松弛迭代法(SOR法)。由Jacobi迭代公式x(k 1) Bx(k) g可知，迭代的每一步计算过程，都是用x(k) 的全部分量来计算x(k1)的所有分量，显然在计算第i个分量Xi(k1)时，已经计算出 来的最新分量x,(k 1)，X2(k 1)x 1(

19、k 1)没有被利用，从直观上看，最新计算出来的分量可能比旧的要好些，因此这些最新计算出来的第i+1次近似x(k 1)的分量Xj(k 1)代只需用n个单元储存近似解分替老分量Xj(k)进行计算。这样，在整个计算过程中, 量。选取初始向量x(0)，用迭代矩阵(k 1)D -L -1 Ux(k) + D - L -1Ub产生近似解序列 x(k)，这种方法叫Gauss-Seide迭代法。在Gauss-Seidel迭代法 基础上进行加权平均的一种新的迭代法，即i 1nxi(k 1) (1)x/k)biajXjk1ajXjki 1,2L ,naiij 1j i+1其中 称为松弛因子，当 1时称为低松弛，1

20、时是Gauss-Seide迭代法，当1时称为超松弛法，简称SOR法(Successive Over-Relaxation)松弛法的迭代方 式的矩阵表示为x(k 1)=x(k)+X=1-X(k)+ D 1Lx k1D 1Ux kD 1b因为ID 1L1，故1D1 -1L存在，从而有DL-1存在，所以X(k 1)= D-L11D Ux(k)+ D1Lb松弛法的迭代矩阵为M1=D- L1D U本问题使用Gauss-Seide迭代法与SOR法。3. 算法本题在程序实现上使用的是矩阵迭代法，所以Gauss-Seide迭代法与SOR迭代法在程序算法上只有迭代公式不同，下面只介绍Gauss-Seide迭代法的算法。(1) 输入 A=4,3,0;3,4,-1;0,-1,4;b=24;30;-24;x0=(1,1,1); 维数 3，以及迭代次数 7。(2) 计算迭代矩阵B,以及g,进而计算出初次迭代X。若i7，输出结果X，停机，否则转3。4. 程序(1)Gauss-Seide迭代法(2) SOR法(3 =1.25clearclci=1;A=4,3,0;3,4,-1;0,-1,4; D=diag(diag(A);L=D-t