《计算方法的程序设计》由会员分享，可在线阅读，更多相关《计算方法的程序设计（42页珍藏版）》请在金锄头文库上搜索。
1、计 算 方 法 （B）上机实习报告姓 名： 李晋学 号： 3107162006 学 院： 电 气 学 院 考核方式： 考 试以下程序都采用win-Tc进行编译的3-564-2-3811-9151-922-175-1611A=-11327-1-2431-721129-811-1-4-172-127-19一、用列主元GAUSS消去法解线性方程组AX=b，其中:1.1算法组织GAUSS 消去及回代计算b=1122995825用途：解一个n x n的线性方程组输入：方程个数n，系数矩阵A,右端向量b 输出：方程组的解 x , x , , x ，或者报告方程组没有唯一解12 n算法 GAUSSPP(A,
2、b,n,x)本算法按列主元素法求n阶线行方程组Ax=b的解，数据组织同GAUSS1 For k=1,2, .,n-11.1找满足丨aukk | =max | aik丨的下标uk1.2 If aukk=0 then 输出错误信息;stop1.3 For j=1,2,. ,n1.3.1 a . o a ,.kj ukj1.3.2 0 _ o B _kuk1.4 For i=k+1,k+2,. ,n1.4.1 aik/akk=aik1.4.2 For j=k+1,k+2,. ,n1.4.2.1. a -a a =aij ik kj ij1.4.3 0 i-aik0 k=0 i2.回代过程2.10 /
3、a =xn nn n2.2 For k=n-1, n-2,. ,1221 (0 k- akjxj )/akk=Xkj=k+11.2 程序Gauss 列主元#include#include#includevoid GAUSSPP(double a77,double b,int n,double x)int i,j,k,m;double amax,tempt,s=0;for(k=0;kn;k+)amax=akk;m=k;for(i=k+1;ifabs(amax)amax=aik;m=i;if(amk=0) printf(error!);break;if(m!=k) for(j=0;jn;j+)te
4、mpt=akj;akj=amj;amj=tempt;tempt=bk;bk=bm;bm=tempt;for(i=k+1;in;i+)aik=aik/akk;for(j=k+1;j=0;k-) s=bk;for(j=k+1;jn;j+)s=s-akj*xj;xk=s/akk; /*The result is out*/main() double a77=3,-5,6,4,-2,-3,8,1,1,-9,15,1,-9,2,2,-1,7,5,-1,6,11,-1,1,3,2,7,-1,-2,4,3,1,-7,2,1,1,2,9,-8,11,-1,-4,-1,7,2,-1,2,7,-1,9;doubl
5、e b=11,2,29,9,5,8,25;double x7=0,0,0,0,0,0,0;int i;clrscr();GAUSSPP(a,b,7,x);printf(The result is X= n); for(i=0;i7;i+)printf(x%d=%fn,i+1,xi); getch();1.3计算结果及结果分析GAUSSPP* The result is X=*x1=1.000000x2=1.000000x3=1.000000x4=1.000000x5=1.000000x6=1.000000 x7=1.000000此题中高斯消去法能够顺利进行是因为每一步都满足条件akk(k-i)
6、工0。所 编程序均米用双精度浮点运算,运算精度较高，而且GAUSS列主兀GAUSS消去法 是稳定的算法，矩阵的条件数cond(A)都较小，所以问题是良态的，而且系数矩阵的 扰动A为零。但是实际上，实数在计算机中，以二进制浮点数存在。浮点数是 离散的有限数集。十进制实数在输入计算机时，有可能因为进制转换和二进制浮 点精度有限的原因产生误差，例如0.1， 1/3这两个十进制实数。在运算中，也会 积累舍入误差。上面显示结果没有误差，是因为计算程序中使用了 64 位的 double 变量，所以输入和计算过程中，最终误差远小于1.0e - 8。二、矩阵LDLT分解与Cholesky分解：求矩阵A=(a
7、)的LDLt分解,及其Cholesky分解，其中ij 20*20”i = ja .= ij min(i, j) i 丰 j2.1 算法组织LDLT 分解用途：把一个对称正定矩阵分解为LDLt的形式，其中L是一个下三角矩阵，其 对角线兀素全部为 1， D 是一个兀素全部为正数的对角线矩阵。输入：矩阵的行数n，对称正定矩阵A = (a ), 1 i n， 1 j nij输出：矩阵L中对角线以下的元素l 1 j Vi 1 i n ，D的元素d， 1 i rpp kp p1.2 a -艺 a r = akk kp p kkp=11.3 If d=0 thenk1.3.1 输出错误信息；stopelse
8、1.3.2 For I=k+1,k+2, ,n1.3.2.1 (a -E a r )/ a = aik kp p kk ikp=11.3.2.2 (a -Em a ) /a = aki p pi kk ik p=1Cholesky 分解用途：把一个对称正定矩阵分解为GGT其中G是一个下三角矩阵 输入：矩阵的行数n，对称正定矩阵A = C ), 1 i n， 1 j nij输出：矩阵G中对角线及其以下的元素l 1 j i 1 i aijip jp jjijp=11.2 (a -Sa a )1/2 = a )iiip ipiip=12.2 程序( 1 ) LDLT#include #include
9、 #define N 20void ggt(double aNN,int n)int i,j,k;double m; for(i=0;in;i+)for(j=0;ji;j+)m=0;for(k=0;kj;k+)m+=aik*ajk; aij=(aij-m)/ajj;m=0;for(k=0;ki;k+)m+=aik*aik;aii=sqrt(aii-m);main()int i,j,n=N;double aNN,lNN,dNN,gNN; for(i=0;in;i+)for(j=0;jn;j+)if(i=j) aij=lij=gij=i+1; else aij=lij=gij=(ij?i:j)+1
10、; ggt(g,n); for(i=0;in;i+) for(j=0;jn;j+) if(ij) gij=0;G=);printf(nnn for(i=0;in;i+) printf(n); for(j=0;jn;j+) printf( %3.0f,gij);getch();(2) Cholesky #include #include #define N 20void ggt(double aNN,int n)int i,j,k; double m; for(i=0;in;i+)for(j=0;ji;j+)m=0;for(k=0;kj;k+) m+=aik*ajk;aij=(aij-m)/aj
11、j; m=0; for(k=0;ki;k+) m+=aik*aik; aii=sqrt(aii-m);main()int i,j,n=N;double aNN,lNN,dNN,gNN; for(i=0;in;i+)for(j=0;jn;j+)if(i=j) aij=lij=gij=i+1; else aij=lij=gij=(ij?i:j)+1; ggt(g,n); for(i=0;in;i+) for(j=0;jn;j+) if(ij) gij=0;G=);printf(nnn for(i=0;in;i+) printf(n); for(j=0;jn;j+) printf( %3.0f,gij);getch();2.3 计算结果及结果分析:111111111d苗苗苗 !1
