计算方法实验

1、数值计算方法专业:姓名:姓名:实验一1. 题目(课本46页习题二第1题)用列主元素法解下列方程，结果保留4位小数。2x16x24x34x14x25x336x1x218x322. 模型建立(理论分析)用方程组的增广矩阵a11a12a1nb1A,ba21a22a2nb2an1an2annbn表示，并在增广矩阵上进行计算。首先在矩阵的第一列中选取绝对值最大的元，通过行互换，将最大的元所在行换在第一行的位置，得到增广矩阵A1,b1，然后进行消元，得到矩阵22A ,b11a11a122 a221 a1n2 a2nb11b22.an22ann然后在每次从消元之前，都通过行变换将绝对值最大的元所在的行换到相

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;1,4,-5;6,-1,18;b=4;3;2;B=Ab;RA=rank(A)

3、;RB=rank(B);zhica=RB-RA;ifzhica>0,disp('请注意：因为RA=RB,所以此方程组无解)returnendifRA=RBifRA=3disp('请注意：因为RA=RB=3,所以此方程组有唯一解)X=zeros(3,1);C=zeros(1,4);forp=1:2Y,j=max(abs(B(p:3,p);C=B(p,:);B(p,:)=B(j+p-1,:);B(j+p-1,:)=C;fork=p+1:3m=B(k,p)/B(p,p);B(k,p:4)=B(k,p:4)-m*B(p,p:4);endendb=B(1:3,4);A=B(1:3,

4、1:3);X(3)=b(3)/A(3,3);forq=2:-1:1X(q)=(b(q)-sum(A(q,q+1:3)*X(q+1:3)/A(q,q);endXelsedisp（'请注意：因为RA=RB<3,所以此方程组有无穷多解）endend5. 运行结果,实验二1. 题目（课本46页习题二第3题）用紧凑格式解下列方程组，并写出L，U矩阵。1234x1214916x210182764x34411681256x41902. 模型建立（理论分析）已知用矩阵形式表示Gauss消去法的消元过程是对方程组增广矩阵A,b进行一系列初等变换，将系数矩阵A左乘上三角形矩阵的过程，也等价于用一串初

5、等矩阵去做乘增广矩阵，因此消元过程可以通过矩阵运算来表示。例如，第一次消元等价于用初等矩阵1000l211 00L1l31 0 10左乘矩阵 A1 ， b1 A， b.ln1001其中li1ai11/a111i2,3,n，即A2，b2L1A1，b1。一般地，第k次消元等价于用初等矩阵101lk1k1左乘矩阵 A k ， b k ，其中 liklnk0aikk/akkkik1,n，经过n1次消元后得到An，bnLn1Ln2L1A1，b1，因为Lk11lk1k1lnk0111LL11L21Ln111l211l31l321A，bLAn，Lbln1ln2ln3lnn1消元过程实际上是把系数矩阵A分解成

6、单位下三角与上三角矩阵的乘积的过程。ALU其中L为单位上三角矩阵，u1121311l32u12u22u13u23u1nu2nln1ln2ln31lnn1unn这种分解称为杜利特尔(Doolittle)分解，也称为LU分解。A为n阶方阵，若A的顺序主子式Ai一的Doolittle分解。由矩阵的乘法法则，得由矩阵的乘法法则,i1,2,n得lj和Uj1均不为零，则矩阵A存在唯的公式u1ja1jj1,2,nuuijaaijlijaaiji1likukjk1j1likukjk1i2,/ujjj计算过程应按第1行，第1列，第2行，第再由Lyb和Uxy计算方程组的解。本题要求用紧凑格式解下列方程组，行Doo

7、little分解。,n,ji,1,2,n2列，,n1,ij1,n的顺序进行。得到L，U矩阵并写出L，U矩阵。在程序实现上就是进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的储存空间。u1ja1jj1,2,ni1(2)计算uijaijlikukji2,n,ji,n。k1j1lijaijlikukj/ujjj1,2,n1,ij1,nk1(3) Lyb和Uxy计算方程组的解，停机。4. 程序LU矩阵分解clearclcA=1,2,3,4;1,4,9,16;1,8,27,64;1,16,81,256;b=2

8、;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);fork=2:nu(k,k:n)=A(k,k:n)-l(k,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);endulY=zeros(n,1);Y(1)=b(1);fori=2:nY(i)=b(i)-l(i,1:i-1)*Y(1:i-1);endX=zeros(n,1);ifdet(u)=0;X=0;elseX(n)=Y(

9、n)/u(n,n);fori=n-1:-1:1X(i)=(Y(i)-u(i,i+1:n)*X(i+1:n)/u(i,i);endendX5. 运行结果4-MATLAB7.6.0(R2008s)匕AleEditDebugParallelDesktop'AindowHelp图句TM国明-2014200093匚|£,Command ' Vindow' Vorkspate.123402120062400024iiao13101761OVR .：Start实验三1. 题目（课本47页习题二第6题）分别用平方根法和改进平方根法求解方程组1 213x112 505x2210

10、141x3163 5115x482. 模型建立（理论分析）依据定理，如果A为n阶对称正定矩阵，则存在惟一的非奇异下三角阵L使得ALLT且L的对角元素皆为正数。l11l 21 l 22Al11l 21l 22l n1l n2l n1ln2l nnlnn其中 lii 0 i 1,2 n 。由矩阵乘法及lik 0 （当 ik 时），得naij l ik l jkk1j1lik l jk l jj l ij k1矩阵的这种分法称为Cholesky分解法，即平方根分解法。比较A与LLT的相应元素，可得1/2k12lkkakklkjk1,2,nj1k1likaiklijlkj/lkkik1,nj10规定0

11、。计算顺序是按列进行的，即j1l11 li1 i 2,3, ,nl22li2 i 3, n当矩阵A完成Cholesky分解后，求解对称正定方程组Axb的解计算公式k1ykbklkjyj,k1,2,nj1nXk(ykljkXj)lkk,kn,n1,1jk1当A为对称矩阵，且A的所有顺序主子式均不为零,A可唯一分解为如下的.形式ALU为了利用A的对称性，将U再分解为:U12u11U11U22u nnU1 nu11un-1,nu n-1,n-11DU。其中D为对角阵，U。为单位上三角阵，于是ALULDU0由分解的唯一性即得UoLT,代入得到对称矩阵A的分解式ALDLTo这是一种改进的平方根法（LDL

12、T法）其中L为单位下三角阵，Ddiagd1,dn为对角阵。由比较法可以导出LDLT分解的计算公式：对于k1,2,nk1dkakklkjdjj1k1likaikIjdjlkj/dkik1,nj1计算顺序如下：d1Ii1i2,3,nd2Ii2i3,nda此时式中有许多计算式重复的，改进辅助量Uiklikdkk1,2,n,ik1,nk1dkakkUkjlkjj1k1ukjaikuijlkjj1likUkj/dkLDLT分解后，解对称正定方程组Axb可分两步进行：先解方程组Lyb,再由LTxD1y求x。本问题使用平方根法与改进的平方根法两种方法解方程组。3.算法两种解法算法应相似，只介绍平方根法。(1

13、)输入A=1,2,1,-3;2,5,0,-5;1,0,14,1;-3,-5,1,15;b=1;2;16;8;.lkk(2)计算likakkk1lk2jj11/2k1,2,k1aiklijlkj/lkkikj1,n，得到L矩阵1,n(3)求解先解方程组Lyb,再由LTxy求x。4. 程序(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);forj=2:n-1L(j,

14、j)=sqrt(A(j,j)-sum(L(j,1:j-1).A2);fori=j+1:nL(i,j)=(A(i,j)-sum(L(i,1:j-1).*L(j,1:j-1)/L(j,j);endendL(n,n)=sqrt(A(n,n)-sum(L(n,1:n-1).A2);Lx=zeros(n,1);y=zeros(n,1);y(1)=b(1)/L(1,1);fori=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);fori=n-1:-1:1x(i)=(y(i)-sum(L(i,i+

15、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;8;L=zeros(4,4);U=zeros(4,4);D=zeros(4,4);fori=1:4L(i,i)=1;endD(1,1)=A(1,1);fori=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);fori=3:4U(i,2)=A(i,2)-U(i,1)*L(2,1);L(i,2)=U(i,2

16、)/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,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);fori=2:4y(i)=(b(i)-sum(L(i,1:i-1)'.*y(1:i-1)/L(i,i);endyL=L'y=inv(D)*y;x(4

17、)=y(4)/L(4,4);fori=3:-1:1x(i)=(y(i)-sum(L(i,i+1:4)'.*x(i+1)/L(i,i);endx5. 运行结果(1)平方根法求解方程组程序运行结果九MATLAB7.6.0(R200Sa1早WI,实验四1 .题目(课本72页习题三第4题)分别用Gauss-Seide迭代法与SOR法(=1.25求解方程组430xi2434-1X2300-14x324取x0111T,迭代7次，并比较它们的计算结果。2 .模型建立(理论分析)迭代法主要解决大型稀疏矩阵方程组，该方法包括Jacobi迭代法，Gauss-Seidel迭代法，超松弛迭代法(SOR法)。由

18、Jacobi迭代公式x(k1)Bx(k)g可知，迭代的每一步计算过程，者B是用x(k)的全部分量来计算x(k1)的所有分量，显然在计算第i个分量xjk1)时，已经计算出来的最新分量x,(k1),x2(k1)x:1)没有被利用，从直观上看，最新计算出来的分i+1次近似x(k 1)的分量xj(k1)代只需用n个单元储存近似解分量可能比旧的要好些，因此这些最新计算出来的第替老分量xj(k)进行计算。这样，在整个计算过程中,量。选取初始向量x(0),用迭代矩阵产生近似解序列x(k)x(k 1) D - L -1Ux(k)-1+ D-L Ub,这种方法叫Gauss-Seide迭代法。在Gauss-Sei

19、del迭代法基础上进行加权平均的一种新的迭代法，即i 1xi(k1) (1)xi(k) - biaj xjk 1aiij 1其中 称为松弛因子，当 1时称为低松弛，nkaijxji 1,2L ,nj i+11时是Gauss-Seidel迭代法，当1时称为超松弛法，简称SOR法(Successive Over-Relaxation)松弛法的迭代方 式的矩阵表示为x(k 1)=x(k) +1-x(k)1, k 11, k1,x + D Lx D Ux D b因为I D1L 1,故I D1L-1存在，从而有D L-1存在，所以x(k 1)= D-松弛法的迭代矩阵为M = D- LD U x(k)+

20、D L 1 b1本问题使用Gauss-Seide迭代法与SOR法3 .算法本题在程序实现上使用的是矩阵迭代法，所以Gauss-Seidel迭代法与SOR迭代法在程序算法上只有迭代公式不同,下面只介绍Gauss-Seide迭代法的算法。(1)输入A=4,3,0;3,4,-1;0,-1,4;b=24;30;-24;x0=(1,1,1);维数3，以及迭代次数7。计算迭代矩阵B,以及g,进而计算出初次迭代x。(3)若i<=7，令x0=x，再次迭代，并记循环次数i=i+1。(4)若i>7，输出结果X，停机，否则转3。4. 程序(1) Gauss-Seidel迭代法clearclci=1;A=

21、4,3,0;3,4,-1;0,-1,4;D=diag(diag(A);L=D-tril(A);U=D-triu(A);b=24;30;-24;x0=ones(3,1);X=zeros(7,3);B=inv(D-L)*U;g=inv(D-L)*b;x=B*x0+g;while i<=7x0=x;X(i,:)=x0'x=B*x0+g;i=i+1;endX(2) SOR法(=1.25clearclci=1;w=1.25;A=4,3,0;3,4,-1;0,-1,4;D=diag(diag(A);L=D-tril(A);U=D-triu(A);b=24;30;-24;x0=ones(3,1

22、);X=zeros(7,3);M=inv(D-w*L)*(1-w)*D+w*U);g=inv(D-w*L)*w*b;x=M*x0+g;while i<=7x0=x;X(i,:)=x0'x=M*x0+g;i=i+1;endX5. 运行结果(1) Gauss-Seidel迭代法运行结果4MATLAB7,6,0(R2OO8a)1口gFileEditDebugP5rjllelDeslktopVindowHelpda总弓d*i才苴QC：UsetsAdministrato|.AjCommandIVind5nWorkspace5.25003.8125-5.04693.14063,8828-B.

23、02933.08793.92S8-5.01333.05493,9542-5.011413.03433.5714-5,007213.0215上9821-i00453.01343.95S8-5.0028»4 StartOVR(2) SOR法(=1.25运行结果,.实验五1 .题目（课本139页习题五第3题）已知函数表如下:X1346fXi-75814试求其3阶Lagrange插值多项式，并以此计算f2的近似值。2 .模型建立（理论分析）当有n+1个节点时，构造n次多项式li xli Xji0,1,2,n,它们满足ii,.然后以相对应点的函数值为系数线性组合,lixi0,1,2,n的表达式

24、。此时的lili1,故它必是以下多项式既得所要求的插值多项式。下面推导X有n个根Xjj0,1,2,n,ji,切liXxXoXXonXxX1xXi1XiX1XXi1X%XiXnj0XiXji0,1,nXj这种函数称为Lagrange插值基函数。利用它们得出nnXyJiXi0nVi0j0XiXj此式称为n次Lagrange插值多项式。常写为nLnxyJii0本题要求计算3阶Lagrange插值多项式,插值多项式。并用于计算非节点处的函数近似值。即4个节点的n3次的Lagrange3.算法(1)输入x0=1,3,4,6;y0=-7,5,8,14;并设函数量x。(2)对于linXxjLi0,1,n。j

25、0XiXjnLyJi。i0(4)输出Lfx,停机。4 .程序Lagrange插值程序段clearclcx0=1,3,4,6;y0=-7,5,8,14;symsx;n=size(x0,2);sum=0;fori=1:n;l=1;forj=1:n;ifj=i;l=l*(x-x0(j)/(x0(i)-x0(j);endendsum=sum+l*y0(i);endL=simplify(sum)5 .运行结果MATLAB7.6.0(RZOOGa)|口至尸二i实验六1 .题目（课本140页习题五第8题）已知函数表如下:Xii.6i5i.634i.702i.828fXi2.4i4502.462592.652

26、7i3.03035试求其3阶Newton插值多项式，并以此计算fi.682的近似值。2 .模型建立（理论分析）线性插值公式可表示为iXfX0XX0fX0,Xi称为一次Newton插值多项式。一般地，由各阶差商的定义，依次可得XfX0XX0fX,X0X,X0f X0,XiX Xif x, X0,XiX,X0,XiX0, Xi,X2X X2 f X, X0, Xi,X2X,X0,Xn if X0,Xi,XnX Xn f X,X0,Xn将以上各式分别乘以i, XX0 , X X0 X Xi , X X0X XiX Xn i 然后两边再消去相等的部分，即得n次Newton插值多项式。Nn Xf X0X X0 f X0,XiX X0 X XiX Xn iX X0 X X1f f X0,Xi,