计算方法数值分析C语言源程序.docx

《计算方法数值分析C语言源程序.docx》由会员分享，可在线阅读，更多相关《计算方法数值分析C语言源程序.docx（109页珍藏版）》请在三一文库上搜索。

1、-精品文档!值得拥有！-第1章 线性方程组的直接算法求解线性方程组的直接算法是基于矩阵分解的算法。常见的矩阵分解有两种：矩阵的三角分解矩阵的三角就是把一个矩阵分解成两个三角形矩阵的乘积。比如：简单的三角分解：，这里，是单位下三角矩阵，是上三角矩阵。列主元三角分解：，这里，是初等置换阵，是单位下三角矩阵且各元素的模不超过1，是上三角矩阵。全主元三角分解：，这里，是初等置换阵，是单位下三角矩阵且元素的模不超过，是上三角矩阵。矩阵的正交三角分解正交化三角化就是把一个矩阵分解成一个正交矩阵和一个上三角矩阵的乘积即，这里 是正交矩阵，是上三角矩阵1.1 矩阵的三角分解1.1.1 功能把实矩阵分解成单位下

2、三角形矩阵和上三角形矩阵的乘积即该算法适用于各阶顺序主子式不等于的矩阵1.1.2 算法概述所谓三角分解就是把阶方阵作如下分解其中是单位下三角矩阵，上三角矩阵。当时，构造Gauss变换则以此类推，只要对角线上的元素，便可以一直这样做下去。直到将其化为上三角矩阵为止。即有其中，且从而有而且综上所述，我们可以将两个矩阵因子继续存储在原矩阵的存储空间上的主对角线上的不予存储算法5.3（计算三角分解：Gauss消去法）1.1.3 算法程序1.1.3.1.1 参数说明*a n阶矩阵； n 矩阵的阶；1.1.3.1.2 程序bool GaussLU(double *a,int n)/n阶矩阵的LU分解for

3、(int k=0;kn-1;k+)if(akk=0)return false;for(int i=k+1;in;i+)aik/=akk;for(int j=k+1;jn;j+)aij-=aik*akj;return true;1.1.4 例题求矩阵1 2 3 41 4 9 161 8 27 64的三角分解，结果如下：1 2 3 4 1 2 6 12 1 3 6 24 1 7 6 24 1.2 列主元三角分解1.2.1 功能用矩阵的列主元三角分解，分解矩阵：，这里，是初等置换阵，是单位下三角矩阵且各元素的模不超过1，是上三角矩阵。1.2.2 算法概述列主元三角分解法和普通三角分解法基本上类似，所

4、不同的是在构造Gauss变换前，先在对应列中选择绝对值最大的元素（称为列主元），然后实施初等行交换将该元素调整到矩阵对角线上。例如第步变换叙述如下：选主元：确定使；行交换：将矩阵的第行和第行上的元素互换位置。即实施Gauss变换：通过初行变换，将列主对角线以下的元素消为零即算法5.2（计算三角分解：列主元Gauss消去法）1.2.3 程序设计1.2.3.1 参数说明*a：n阶矩阵； *p：记忆分解过程中进行的行交换； n：矩阵的阶1.2.3.2 程序int RowGaussLU(double *a,int *p,int n)/n阶矩阵的列主元LU分解int i,j, k;double max;

5、for(i=0;in;i+) pi=i;for(i=0;in-1;i+)max=aii;pi=i;for(j=i+1;jfabs(max)max=aji;pi=j;if(max=0) return i-1;else if(pi!=i)for(j=0;jn;j+)swap(aij,apij);for(j=i+1;jn;j+)aji/=aii;for(j=i+1;jn;j+)for(k=i+1;kn;k+)ajk-=aik*aji;if(an-1n-1=0) return n-1;elsereturn n;1.2.4 例题用列主元高期分解矩阵结果，P=(0,3,2,3)1.3 全主元三角分解1.3

6、.1 功能将矩阵进行如下三角分解，其中：是初等置换阵，单位下三角阵（各元素绝对值不超过），上三角阵1.3.2 算法概述全主元三角分解法和普通三角分解法基本上类似，所不同的是在构造Gauss变换前，先在对应列中选择绝对值最大的元素（称为列主元），然后实施初等行交换将该元素调整到矩阵对角线上。例如第步变换叙述如下：选主元：确定，使；行列交换：将矩阵的第行和第行上的元素互换位置，再将第列和第列上的元素互换。即；实施Gauss变换：通过初行变换，将列主对角线以下的元素消为零即算法5.2（计算三角分解：列主元Gauss消去法）1.3.3 算法程序1.3.3.1 参数说明*a: 待分解的目标矩阵，也是结果

7、存储的实体 *P: 保存行主元行交换次序*q：保存主元列交换次序； n：矩阵的阶1.3.3.2 程序int AllGaussLU(double *a,int *p,int *q,int n)/n阶矩阵的全主元LU分解int i,j,k;double max;for(i=0;in;i+)pi=qi=i;for(i=0;in-1;i+)max=aii;pi=i;qi=i;for(j=i;jn;j+)for(k=i;kfabs(max)max=ajk;pi=j;qi=k;if(max=0) return i-1;if(pi!=i)for(j=0;jn;j+)swap(aij,apij);if(qi!

8、=i)for(j=0;jn;j+)swap(aji,ajqi);for(j=i+1;jn;j+)aji/=aii;for(j=i+1;jn;j+)for(k=i+1;kn;k+)ajk-=aik*aji;if(an-1n-1=0) return n-1;elsereturn n;1.3.4 例题用全主元三角分解矩阵1 2 3 41 4 9 161 8 29 641 16 81 256系数矩阵分解为256 81 1 16 0.25 6.75 0.75 4 0.015625 0.256944 0.791667 0.722222 0.0625 0.583333 0.631579 0.210526 全

9、主元选取信息：行信息P：3 2 3 3 列信息Q：3 2 3 3 1.4 对称正定Cholesky分解1.4.1 功能将对称正定矩阵能作以下分解：，其中为下三角形矩阵，为的转置。1.4.2 算法概述，以列为序从左到右，用待定系数法，则的元素（可由下列公式确定：另一方面，则的元素（也可由Gauss消去法的思想逐列得到，算法如下：1.4.3 算法程序1.4.3.1 参数说明*a：对称正定矩阵，算法结束后下三角部分保存Cholesky因子；n：矩阵的阶1.4.3.2 程序void CholeskyLL(double *a,int n)/n阶对称正定矩阵的LL分解int i,j,k;for(k=0;k

10、n;k+) akk=sqrt(akk);for(i=k+1;in;i+)aik/=akk;for(j=k+1;jn;j+)for(i=j;in;i+)aij-=aik*ajk;1.4.4 例题对称正定矩阵：，Cholesky因子1.5 改进的Cholesky三角分解1.5.1 功能该算法将阶实对称正定矩阵分解为：，其中为单位下三角形矩阵，是对角矩阵，为的转置。1.5.2 算法概述设 ，以列为序从左到右，用待定系数法，则的元素（以及对角矩阵的元素可由下列公式确定：1.5.3 算法程序1.5.3.1 参数说明*a：对称正定矩阵，算法结束后下三角部分保存Cholesky因子；n：矩阵的阶1.5.3.

11、2 程序void CholeskyLDL(double *a,int n)/n阶对称正定矩阵的LDL分解double *v=new doublen;for(int j=0;jn;j+)for(int i=0;ij;i+)vi=aji*aii;for(i=0;ij;i+)ajj-=aji*vi;for(i=j+1;in;i+) for(int k=0;kj;k+)aij-=aik*vk;aij/=ajj;delete v;1.5.4 例题对称正定矩阵：，分解结果：1.6 下三角形方程组1.6.1 功能该算法用来求解下三角形方程组。它又称为前代法。1.6.2 算法概述求解上三角形方程组的前代法，我

12、们通常采用顺序消元法实现。具体地，从其增广矩阵依次沿下述步骤进行：第1步：消去后个方程中的，即其中，。第2步：再消去后个方程中的，即其中，。直到第步：消去第个方程中的，即其中，第步：计算。到此我们得到与原上三角方程组同解的对角形方程组，其增广矩阵为其解是显然的，即。1.6.3 算法分析算法5_1（前代法：解上三角方程组）当系数矩阵是单位下三角矩阵时，算法简化为beginfor k=1:n-1for i=k+1:nb(i)=b(i)-a(i,k)b(k)endendend1.6.4 程序设计/前代法void Lower(double *a,double *b,int n) for(int i=0

13、;in-1;i+)bi/=aii;for(int j=i+1;jn;j+)bj-=bi*aji;bn-1/=an-1n-1;/单位下三角方程组前代法void UnitLow(double *a,double *b,int n) for(int i=0;in;i+)for(int j=i+1;j0;i-)bi/=aii;for(int j=0;j=0;i-)for(int j=0;ji;j+)bj-=bi*aji;1.8 高斯消消法1.8.1 功能该算法用于求解简单的线性方程组其中系数矩阵须满足顺序主子式不等于零。1.8.2 算法概述1.8.2.1 算法基于系数矩阵的三角分解。算法计算步骤如下：

14、将作三角分解：；求解下三角方程组：；求解上三角方程组：1.8.3 算法程序1.8.3.1 参数说明*a：n阶矩阵； *b：常数项； n：矩阵的阶1.8.3.2 程序bool Gauss(double*a,double*b,int n)If(GaussLU(a,n)UnitLow(a,b,n);Uper(a,b,n);Return true;elsereturn false;1.8.4 例题用高期消去法解线性方程组系数矩阵的三角分解，结果如下：1 2 3 4 1 2 6 12 1 3 6 24 1 7 6 24 下三角方程组的解：(2，8，18，34)上三角方程组的解（即原方程组的解）：（-1，

15、1，-1，1）1.9 解三对角型方程组的追赶法1.9.1 功能本过程用以求解三对角型方程组：（1）1.9.2 方法概要首先对A作下列分解：（2）其中矩阵和的形式如下：，而和的元素由下列关系式确定：（）进而解方程组（）就等价于解方程组：（4）（5）其计算公式如下：（6）（7）（6）式为“追”的过程，（7）式为“赶”的过程，故称该法为追赶法。1.9.3 程序设计1.9.3.1 参数说明*a 存储三对角矩阵。a0: 下次对角线；a1：主对角线；a2：上次对角线。*b 右端向量； 方程组的阶1.9.3.2 C程序void Pursuit(double *a,double *b,int n)/追赶法do

16、uble w;int k;w=a10;b0/=a10;for(k=0;k=0;k-)bk=bk-a1k*bk+1;1.9.4 例题求解三对角型方程组：计算结果。1.10 列主元高斯消去法1.10.1 功能用列主元高斯消去法求解线性方程组1.10.2 算法概述算法基于列主元三角分解 计算步骤如下：用列主元三角分解法（算法5_2）对系数矩阵进行列主元三角分解：即对向量进行行调整：；求解单位下三角方程组：；求解上三角方程组：。1.10.3 程序设计1.10.3.1 参数说明*a：n阶矩阵； *b：常数项； n：矩阵的阶1.10.3.2 程序bool RowGauss(double *a,double

17、 *b,int n)/列主元高斯消去法int *p=new intn;if(RowGaussLU(a,p,n)n)return false;for(int i=0;in-1;i+)if(pi!=i)double t=bi;bi=bpi;bpi=t;UnitLow(a,b,n);Uper(a,b,n);return true;1.10.4 例题用列主元高期消去法解线性方程组系数矩阵分解为：，其中，；的解为：=（2，188，-38.5714，2.18182）；的解，即原方程组的解：（-1，1，-1，1）。1.11 全主元高斯消去法1.11.1 功能用全主元Gauss消去法求解线性方程组1.11.2

18、 算法概述算法基于全主元三角分解计算步骤如下：用列主元三角分解法（算法5_2）对系数矩阵进行列主元三角分解：即对向量进行行调整：；求解单位下三角方程组：；求解上三角方程组：；调整解向量：。1.11.3 算法程序1.11.3.1 参数说明*a：系数矩阵并保存全主元三角分解结果；*b：常数向量； n：方程的阶1.11.3.2 程序bool AllGauss(double *a,double *b,int n)/全主元高斯消去法int *p=new intn;int *q=new intn;if(AllGaussLU(a,p,q,n)n)return false;for(int i=0;i=0;i-

19、)if(qi!=i)double t=bi;bi=bqi;bqi=t;return true;1.11.4 例题用全主元高期消去法解线性方程组系数矩阵分解为256 81 1 16 0.25 6.75 0.75 4 0.015625 0.256944 0.791667 0.722222 0.0625 0.583333 0.631579 0.210526 全主元选取信息：行信息P：3 2 3 3 列信息Q：3 2 3 3 向量调整后：P的解：的解：列调整：1.12 对称正定方程组的平方根法1.12.1 功能平方根法用于求解对称正定线性方程组其中为阶实对称正定矩阵。1.12.2 算法概述解正定方程组

20、不需要选主元，这是由对称正定矩阵的性质所决定的。求解正定方程组的具体计划如下：将系数矩阵A作Cholesky分解：调用Cholesky LL(a,n);调用求解下三角方程组的前代法，求解：调用Lower(a,n)；调用求角上三角方程组的回代法求解：先求A的转置，再调用Uper(a,n)；1.12.3 程序void Cholesky(double *a,double *b,int n)Cholesky LL(a,n);Lower(a,b,n);for(int i=0;in;i+)for(int j=i+1;jn;j+)aij=aji;Uper(a,b,n);1.12.4 例题求解正定方程组计算结

21、果：Cholesky因子：，方程组的解向量为：1.13 改进平方根法1.13.1 功能平方根法用于求解对称正定线性方程组其中为阶实对称正定矩阵。1.13.2 算法概述解正定方程组不需要选主元，这是由对称正定矩阵的性质所决定的。求解正定方程组的具体计划如下：将系数矩阵A作改进Cholesky分解：调用Cholesky LDL(a,n);调用求解下三角方程组的前代法，求解：调用Lower(a,n)；求解对角方程组：调用求角上三角方程组的回代法求解：先求A的转置，再调用Uper(a,n)；1.13.3 算法程序void CholeskyD(double *a,double *b,int n)/正定方

22、程组的平方根算法2CholeskyLDL(a,n);UnitLow(a,b,n);for(int i=0;in;i+)bi/=aii;for(int j=i+1;jn;j+)aij=aji;UnitUp(a,b,n);1.13.4 例题求解正定方程组计算结果：，第2章 最小二乘问题的解法2.1 正则化算法2.1.1 功能通过解法方程组的解求得方程组的最小二乘解。2.1.2 算法概述算法思想基于下面的定理：定理 是的最小二乘解当且仅当。算法1（正则化法：求解的最小二乘解）2.1.3 C程序参数说明：*a：的列满秩矩阵；b：维向量；x：维解向量。void LeastSquares(double *

23、a,double *b,double *x,int m,int n)double *c=new double*n; for(int i=0;in;i+) ci=new doublen;for(i=0;in;i+)for(int j=0;jn;j+)cij=0.0;xi=0.0;for(int k=0;km;k+)cij+=aki*akj;xi+=aki*bk;CholeskyLL(c,n);Lower(c,x,n);for(i=0;in;i+)for(int j=i+1;jn;j+)cij=cji;Uper(c,x,n);for(i=0;in;i+)delete ci;delete c;2.1

24、.4 例题求一个二次多项式使在残向量的2范数最小的意义下拟合下面的-1-0.75-0.500.250.50.751.000.81250.951.001.31251.752.3125计算结果：.2.2 Householder变换2.2.1 功能针对非零向量，构造Householder变换，使2.2.2 算法概述定理 设，则可构造单位向量，使由(3.2.1)定义的Householder变换满足，其中证明 令则定理证明表明，对任意的可按如下的步骤来构造确定：计算；计算构造算法程序时，注意以下几个问题：计算时，为避免溢出，可以考虑先将规范化为的向量。事实上，若，假定Householder变换阵使，那么

25、也就是说，对于向量集合中的任一个非零向量，由Th3.2.2构造出的Householder变换阵都是一样的，或许最多只相差一个符号。在实际计算时，前的符号如何选取最好。即在计算向量的第一个分量时，为了避免两相近数相减，而导致计算结果损失有效数字，可适当做变形处理，即如果时，直接计算：；如果时，事实上，考虑到实际计算的需要，在最后构造Householder变换阵时，常规定H具有如下形式，其中：的第一个分量为1；算法3.2.1 （计算Householder变换）2.2.3 算法程序2.2.3.1 参数说明*x：n维向量 n：向量的维数 返回值：Householder变换中的，*v：Household

26、er变换中向量2.2.3.2 程序double Householder(double*x,double *v,int n)int i;double alph,bet,norm;norm=fabs(x0);for(i=1;in;i+)if(normfabs(xi)norm=fabs(xi);if(norm=0)AfxMessageBox(Fail, because x is zero.n);exit(0);for(i=0;in;i+)xi/=norm;norm=0.0;for(i=1;in;i+)norm+=xi*xi;v0=1.0;for(i=1;in;i+)vi=xi;if(norm=0)b

27、et=0.0;elsealph=sqrt(x0*x0+norm);if(x0=0)v0=x0-alph;elsev0=-norm/(x0+alph);bet=2.0*v0*v0/(norm+v0*v0);alph=v0;for(i=0;in;i+)vi=vi/alph;return bet;2.2.4 例题把向量正交化为与平等且等长的向量的Householder变换为：，且 其中：（1，-0.824343，-0.206086， -0.7213， -0.412171），=0.829128，11.7047.2.3 矩阵的正交三角化2.3.1 功能运用Householder变换做正交约化，将矩阵约化

28、为上梯形矩阵即求得正交矩阵使2.3.2 算法概述定理（QR分解定理）设，则有QR分解：其中是正交矩阵，是具有非负对角元的上三角阵；而且当且非奇异时，上述的分解还是唯一的。用Householder方法计算QR分解与不选主元的Gauss消去法很类似，就是利用Householder变换逐步将约化为上三角矩阵。设，并假定已经计算出Householder变换和使得那么我们现在的任务就是集中精力于第三列杯为“+”的5个元素，确定一个Householder变换使得.令，则有对于一般的矩阵，假定我们已经进行了步，得到了householder变换，使得其中是上三角矩阵。假定第步是：先用算法32.1确定House

29、holder变换使得其中；然后，再计算。令则其中是上三角矩阵。这样，从出发，依次进行次，我们就可将约化为上三角阵。现在记则注意，这样得到的上三角矩阵的对角元均是非负的。下面考虑计算的QR分解的存储问题。当分解完成后，一般来说，就不再需要，便可用来存放和。通常并不是将算出，而是只存放构成它的个Householder矩阵，而对每个，我们只需要保存和即可。注意到有如下形式我们正好可以将存储在的对角元以下位置上。例如，对于的问题，其存储方式如下：综合上面的讨论，可得如下算法。算法3.3.1 (计算QR分解:Householder变换法) 容易算出,这一算法的运算量为2.3.3 算法程序2.3.3.1

30、参数说明A：矩阵，作为输入是待分解的矩阵，作为输出上三角部分存储，下三角部分顺序存储各Householder向量Bet：顺序存储各Householder变换的M和n：指定矩阵的大小2.3.3.2 程序void HouseholderQR(double *a,double*bet,int m,int n)/m*n阶矩阵的QR分解double *x,*v;x=new doublem;v=new doublem;double t;int i,j,k;for(k=0;kn;k+)/初始向量xfor(i=k;im;i+)xi=aik;/调用Householder(),产生向量v和参数betbetk=Ho

31、useholder(x+k,v+k,m-k);/A:=HAfor(j=k;jn;j+)t=0.0;for(i=k;im;i+)t+=vi*aij;t*=betk;for(i=k;im;i+)aij-=t*vi;/存储vfor(i=k+1;im;i+)aik=vi;2.3.4 例题矩阵：1 1 -22 1 01 -1 0-1 2 1QR三角分解结果：(左下三角部分各Gouse变换的v向量)2.64575 0 -1.13389 -1.21525 2.64575 -2.22045e-016 -0.607625 0.911438 1.92725 0.607625 -3.23431 1.25684 各H

32、ouse变换的bet值0.622036 0.162714 0.7753 2.4 求超定线性方程组的正交化算法2.4.1 功能该算法用来求解超定方程组的最小二乘解2.4.2 算法概述设有线性无关的列，并且假定。现将分块为并且令那么由此即知，是最小二乘问题(3.1.3)的解当且仅当是的解。这样一来，最小二乘问题(3.1.3)的解就可很容易从上三角方程组求得。综合上面的讨论，可得正交化方法的基本步骤为：计算的QR分解（3.3.2）；计算；求解上三角方程组.2.4.3 算法程序2.4.3.1 参数说明A：超定方程组系数矩阵 b：超定方程组右端向量 bet：Householder变换顺序值2.4.3.2

33、 程序void LeastSquares(double *a,double *b,double *bet,int m,int n)/最小二乘正交化算法HouseholderQR(a,bet,m,n);double t;for(int j=0;jn;j+)t=bj;for(int i=j+1;im;i+)t+=aij*bi;t*=betj;bj-=t;for(i=j+1;im;i+)bi-=t*aij;Uper(a,b,n);2.4.4 例题用正交化法求方程组的最小二乘解。解结程序处理后得到：增广矩阵正交化后的结果：(左下三角部分各Gouse变换的v向量)11.7047 2.13589 6.15

34、138 -0.824343 3.79973 -5.56321 -0.206086 -0.148331 0.133683 -0.7213 1.36781 2.48331 -0.412171 0.541989 0.16311 各House变换的bet值0.829128 0.627619 最小二乘解：0.79272 -1.46411 第3章 线性方程组的迭代解法本章收集的算法是常见的单步线性定常迭代法即形如其中称为迭代矩阵，是常数项，是初始向量迭代法就是给出了一种构造向量序列的方法，对算法的最基本的要求就是算法的相容性，即迭代方程组（）与原方程组（）是同解的在具体选用算法还需要进一步考虑算法的收敛性

35、，即由算法构造的向量序列是有极限的，而且极限向量恰是方程组的解向量即存在使，且 .单步线性定常迭代法的收敛性充分必要条件是线性方程组的古典迭代算法，通常基于系数矩阵的如下分裂：其中，3.1 雅可比(Jacobi)迭代法3.1.1 功能用雅可比迭代法求解线性方程组的解但算法并不是万能的，其收敛性是有条件的详见教材3.1.2 算法概述雅可比迭代方程为，或迭代矩阵和常数项如下：，雅可比迭代公式，分量形式显然，雅可比迭代要求线性方程组系数矩阵的对角线元素非零并且迭代过程中，与各分量的计算次序无关3.1.3 程序设计3.1.3.1 参数说明A：方程组系数矩阵； b：方程组右端常向量； x：迭代解向量；n

36、：方程个数； Max：最大迭代次数； eps：最大允许误差3.1.3.2 程序double Jacobi(double *A,double *b,double* x,int n,double Max,double eps) double *y=new doublen;double err=1.0; int i,j;for(i=0;i=eps&countMax)count+;err=0.0;for(i=0;in;i+)xi=bi;for(j=0;ji;j+)xi-=Aij*yj;for(j=i+1;jn;j+)xi-=Aij*yj;xi/=Aii;if(errfabs(xi-yi) err=fa

37、bs(xi-yi);for(i=0;in;i+)yi=xi;return count;3.1.4 例题求解线性方程组迭代次数:33近似解向量:x0=0.413793 x1=0.172414 x2=0.862069 误差限: 1e-0063.2 高斯赛德尔(Gauss-Seidel)迭代法3.2.1 功能用高斯赛德尔迭代法求解线性方程组的解但算法并不是万能的，其收敛性是有条件的详见教材3.2.2 算法概述高斯赛德尔迭代方程为，或高斯赛德尔迭代公式，分量形式显然，高斯赛德尔迭代同样要求线性方程组系数矩阵的对角线元素非零但高斯赛德尔迭代与分量的计算次序有关，只能从小到大顺序计算3.2.3 程序设计3.2.3.1 参数说明A：方程组系数矩阵； b：方程组右端常向量； x：迭代解向量；n：方程个数； Max：最大迭代次数； eps：最大允许误差3.2.3.2 程序double GaussSeidel(double *A,double *b,double* x,int n,double Max,double eps) double y,err=1.0; int i,j;for(i=0;i=eps&countMax)count+;err=0.0;for(i=0;in;i