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

《计算方法数值分析C语言源程序1.pdf》由会员分享，可在线阅读，更多相关《计算方法数值分析C语言源程序1.pdf（110页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

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

2、1.1 矩阵的三角分解矩阵的三角分解1.1.11.1.1 功能功能把实矩阵积即该算法适用于各阶顺序主子式不等于的矩阵分解成单位下三角形矩阵和上三角形矩阵的乘1.1.21.1.2 算法概述算法概述所谓三角分解就是把其中是单位下三角矩阵，上三角矩阵。阶方阵作如下分解当时，构造 Gauss 变换则以此类推，只要对角线上的元素样做下去。直到将其化为上三角矩阵为止。即有，便可以一直这其中，且从而有而且综上所述，我们可以将间上两个矩阵因子继续存储在原矩阵的存储空的主对角线上的不予存储算法 5.3（计算三角分解：Gauss 消去法）1.1.31.1.3 算法程序算法程序1.1.3.1.11.1.3.1.1

3、参数说明参数说明*a n 阶矩阵；n 矩阵的阶；1.1.3.1.21.1.3.1.2 程序程序bool GaussLU(double*a,int n)/n 阶矩阵的 LU 分解for(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.41.1.4 例题例题求矩阵1 2 3 41 4 9 161 8 27 64的三角分解，结果如下：1 2 3 41 2 6 121 3 6 241 7 6 241.21.2 列主元三角

4、分解列主元三角分解1.2.11.2.1 功能功能用矩阵的列主元三角分解，分解矩阵：，这里过 1，是初等置换阵，是单位下三角矩阵且各元素的模不超是上三角矩阵。1.2.21.2.2 算法概述算法概述列主元三角分解法和普通三角分解法基本上类似，所不同的是在构造Gauss 变换前，先在对应列中选择绝对值最大的元素（称为列主元），然后实施初等行交换将该元素调整到矩阵对角线上。例如第步变换叙述如下：选主元：确定使行和第；行上的元素互换位置。即行交换：将矩阵的第实施 Gauss 变换：通过初行变换，将列主对角线以下的元素消为零即算法 5.2（计算三角分解：列主元 Gauss 消去法）1.2.31.2.3 程

5、序设计程序设计1.2.3.11.2.3.1 参数说明参数说明*a：n 阶矩阵；*p：记忆分解过程中进行的行交换；n：矩阵的阶1.2.3.21.2.3.2 程序程序int RowGaussLU(double*a,int*p,int n)/n阶矩阵的列主元 LU 分解int i,j,k;double max;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

6、(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.41.2.4 例题例题用列主元高期分解矩阵结果，P=(0,3,2,3)1.31.3 全主元三角分解全主元三角分解1.3.11.3.1 功能功能将矩阵进行如下三角分解，其中：是初等置换阵，单位下三角阵（各元素绝对值不超过），上三角阵1.3.21.3.2 算法概述算法概述全主元三角分解法和普通三角分解法基本上类似，所不同的是在构造Gauss 变换前，先在对应列中选择绝对值最大的元素（称

7、为列主元），然后实施初等行交换将该元素调整到矩阵对角线上。例如第步变换叙述如下：选主元：确定，使；行列交换：将矩阵的第第行和第行上的元素互换位置，再将第列和列上的元素互换。即；实施 Gauss 变换：通过初行变换，将列主对角线以下的元素消为零即算法 5.2（计算三角分解：列主元 Gauss 消去法）1.3.31.3.3 算法程序算法程序1.3.3.11.3.3.1 参数说明参数说明*a:待分解的目标矩阵，也是结果存储的实体*P:保存行主元行交换次序*q：保存主元列交换次序；n：矩阵的阶1.3.3.21.3.3.2 程序程序int AllGaussLU(double*a,int*p,int*q,

8、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!=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+)

9、ajk-=aik*aji;if(an-1n-1=0)return n-1;elsereturn n;1.3.41.3.4 例题例题用全主元三角分解矩阵1 2 3 41 4 9 161 8 29 641 16 81 256系数矩阵分解为256 81 1 160.25 6.75 0.75 40.015625 0.256944 0.791667 0.7222220.0625 0.583333 0.631579 0.210526全主元选取信息：行信息 P：3 2 3 3列信息 Q：3 2 3 31.41.4 对称正定对称正定 CholeskyCholesky 分解分解1.4.11.4.1 功能功能将对

10、称正定矩阵能作以下分解：，其中为下三角形矩阵，为的转置。1.4.21.4.2 算法概述算法概述，以列为序从左到右，用待定系数法，则确定：的元素（可由下列公式另一方面，则算法如下：的元素（也可由 Gauss 消去法的思想逐列得到，1.4.31.4.3 算法程序算法程序1.4.3.11.4.3.1参数说明参数说明*a：对称正定矩阵，算法结束后下三角部分保存 Cholesky 因子；n：矩阵的阶1.4.3.21.4.3.2 程序程序void CholeskyLL(double*a,int n)/n 阶对称正定矩阵的 LL 分解int i,j,k;for(k=0;kn;k+)akk=sqrt(akk)

11、;for(i=k+1;in;i+)aik/=akk;for(j=k+1;jn;j+)for(i=j;in;i+)aij-=aik*ajk;1.4.41.4.4 例题例题对称正定矩阵：，Cholesky 因子1.51.5 改进的改进的 CholeskyCholesky 三角分解三角分解1.5.11.5.1 功能功能该算法将阶实对称正定矩阵分解为：，其中为单位下三角形矩阵，是对角矩阵，为的转置。1.5.21.5.2 算法概述算法概述设，以列为序从左到右，用待定系数法，则的元素的元素（以及对角矩阵可由下列公式确定：1.5.31.5.3 算法程序算法程序1.5.3.11.5.3.1参数说明参数说明*a

12、：对称正定矩阵，算法结束后下三角部分保存 Cholesky 因子；n：矩阵的阶1.5.3.21.5.3.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.41.5.4 例题例题对称正定矩阵：，分解结果：1.61.6

13、下三角形方程组下三角形方程组1.6.11.6.1 功能功能该算法用来求解下三角形方程组。它又称为前代法。1.6.21.6.2 算法概述算法概述求解上三角形方程组的前代法，我们通常采用顺序消元法实现。具体地，从其增广矩阵依次沿下述步骤进行：第 1 步：消去后个方程中的，即其中，。第 2 步：再消去后个方程中的，即其中，直到第，。步：消去第个方程中的，即其中，第步：计算。到此我们得到与原上三角方程组同解的对角形方程组，其增广矩阵为其解是显然的，即。1.6.31.6.3 算法分析算法分析算法 5_1（前代法：解上三角方程组）当系数矩阵是单位下三角矩阵时，算法简化为beginfor k=1:n-1fo

14、r i=k+1:nb(i)=b(i)-a(i,k)b(k)endendend1.6.41.6.4 程序设计程序设计/前代法void Lower(double*a,double*b,int n)for(int i=0;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+

15、)bj-=bi*aji;1.81.8 高斯消消法高斯消消法1.8.11.8.1 功能功能该算法用于求解简单的线性方程组其中系数矩阵须满足顺序主子式不等于零。1.8.21.8.2 算法概述算法概述1.8.2.11.8.2.1 算法基于系数矩阵的三角分解。算法基于系数矩阵的三角分解。算法计算步骤如下：将作三角分解：；求解下三角方程组：；求解上三角方程组：1.8.31.8.3 算法程序算法程序1.8.3.11.8.3.1 参数说明参数说明*a：n 阶矩阵；*b：常数项；n：矩阵的阶1.8.3.21.8.3.2 程序程序bool Gauss(double*a,double*b,int n)If(Gau

16、ssLU(a,n)UnitLow(a,b,n);Uper(a,b,n);Return true;elsereturn false;1.8.41.8.4 例题例题用高期消去法解线性方程组系数矩阵的三角分解，结果如下：1 2 3 41 2 6 121 3 6 241 7 6 24下三角方程组的解：(2，8，18，34)上三角方程组的解（即原方程组的解）：（-1，1，-1，1）1.91.9 解三对角型方程组的追赶法解三对角型方程组的追赶法1.9.11.9.1 功能功能本过程用以求解三对角型方程组：（1）1.9.21.9.2 方法概要方法概要首先对 A 作下列分解：（2）其中矩阵和的形式如下：，而和的

17、元素由下列关系式确定：（）进而解方程组（）就等价于解方程组：（4）（5）其计算公式如下：（6）（7）（6）式为“追”的过程，（7）式为“赶”的过程，故称该法为追赶法。1.9.31.9.3 程序设计程序设计1.9.3.11.9.3.1 参数说明参数说明*a 存储三对角矩阵。a0:上次对角线。*b 右端向量；方程组的阶下次对角线；a1：主对角线；a2：1.9.3.2 C1.9.3.2 C 程序程序void Pursuit(double*a,double*b,int n)/追赶法double w;int k;w=a10;b0/=a10;for(k=0;k=0;k-)bk=bk-a1k*bk+1;1.

18、9.41.9.4 例题例题求解三对角型方程组：计算结果。1.101.10 列主元高斯消去法列主元高斯消去法1.10.11.10.1 功能功能用列主元高斯消去法求解线性方程组1.10.21.10.2 算法概述算法概述算法基于列主元三角分解 计算步骤如下：用列主元三角分解法（算法 5_2）对系数矩阵进行列主元三角分解：即对向量进行行调整：；求解单位下三角方程组：；求解上三角方程组：。1.10.31.10.3 程序设计程序设计1.10.3.11.10.3.1 参数说明参数说明*a：n 阶矩阵；*b：常数项；n：矩阵的阶1.10.3.21.10.3.2程序程序bool RowGauss(double*

19、a,double*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.41.10.4 例题例题用列主元高期消去法解线性方程组系数矩阵分解为：，其中，；的解为：=（2，188，-38.5714，2.18182）；的解，即原方程组的解：（-1，1，-1，1）。1.111.11 全主元高斯消去法全主元高斯消去法1.11.11.

20、11.1 功能功能用全主元 Gauss 消去法求解线性方程组1.11.21.11.2 算法概述算法概述算法基于全主元三角分解计算步骤如下：用列主元三角分解法（算法 5_2）对系数矩阵进行列主元三角分解：即对向量进行行调整：；求解单位下三角方程组：；求解上三角方程组：；调整解向量：。1.11.31.11.3 算法程序算法程序1.11.3.11.11.3.1 参数说明参数说明*a：系数矩阵并保存全主元三角分解结果；*b：常数向量；n：方程的阶1.11.3.21.11.3.2 程序程序bool AllGauss(double*a,double*b,int n)/全主元高斯消去法int*p=new i

21、ntn;int*q=new intn;if(AllGaussLU(a,p,q,n)n)return false;for(int i=0;i=0;i-)if(qi!=i)double t=bi;bi=bqi;bqi=t;return true;1.11.41.11.4 例题例题用全主元高期消去法解线性方程组系数矩阵分解为256 81 1 160.25 6.75 0.75 40.015625 0.256944 0.791667 0.7222220.0625 0.583333 0.631579 0.210526全主元选取信息：行信息 P：3 2 3 3列信息 Q：3 2 3 3向量调整后：P的解：的

22、解：列调整：1.121.12 对称正定方程组的平方根法对称正定方程组的平方根法1.12.11.12.1 功能功能平方根法用于求解对称正定线性方程组其中为阶实对称正定矩阵。1.12.21.12.2 算法概述算法概述解正定方程组不需要选主元，这是由对称正定矩阵的性质所决定的。求解正定方程组的具体计划如下：将系数矩阵 A 作 Cholesky 分解：调用 Cholesky LL(a,n);调用求解下三角方程组的前代法，求解：调用 Lower(a,n)；调用求角上三角方程组的回代法求解Uper(a,n)；：先求 A 的转置，再调用1.12.31.12.3 程序程序void Cholesky(doubl

23、e*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.41.12.4 例题例题求解正定方程组计算结果：Cholesky 因子：，方程组的解向量为：1.131.13 改进平方根法改进平方根法1.13.11.13.1 功能功能平方根法用于求解对称正定线性方程组其中为阶实对称正定矩阵。1.13.21.13.2 算法概述算法概述解正定方程组不需要选主元，这是由对称正定矩阵的性质所决定的。求解正定方程组的具体计划如下：将系数矩阵 A 作

24、改进 Cholesky 分解：调用 Cholesky LDL(a,n);调用求解下三角方程组的前代法，求解：调用 Lower(a,n)；求解对角方程组：调用求角上三角方程组的回代法求解Uper(a,n)；：先求 A 的转置，再调用1.13.31.13.3 算法程序算法程序void CholeskyD(double*a,double*b,int n)/正定方程组的平方根算法 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.41.13.

25、4 例题例题求解正定方程组计算结果：，第第 2 2 章章 最小二乘问题的解法最小二乘问题的解法2.12.1 正则化算法正则化算法2.1.12.1.1 功能功能通过解法方程组的解求得方程组的最小二乘解。2.1.22.1.2 算法概述算法概述算法思想基于下面的定理：定理定理是的最小二乘解当且仅当。算法 1（正则化法：求解的最小二乘解）2.1.3 C2.1.3 C 程序程序参数说明：*a：的列满秩矩阵；b：维向量；x：维解向量。void LeastSquares(double*a,double*b,double*x,int m,int n)double*c=new double*n;for(int

26、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.42.1.4 例题例题求一个二次多项式下面的-0.750.8125使在残向量的 2 范数最小的意义下拟合-1-0.500.251.3

27、1250.50.752.31251.000.951.001.75计算结果：.2.2 Householder2.2 Householder 变换变换2.2.12.2.1 功能功能针对非零向量，构造 Householder 变换，使2.2.22.2.2 算法概述算法概述定理定理 设Householder 变换，则可构造单位向量满足，使由(3.2.1)定义的，其中证明证明 令则定理证明表明，对任意的可按如下的步骤来构造确定：计算；计算构造算法程序时，注意以下几个问题：计算时，为避免溢出，可以考虑先将规范化为的向量。事实上，若，假定 Householder 变换阵使，那么也就是说，对于向量集合中的任一

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

29、older 变换中的，*v：Householder 变换中向量2.2.3.22.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;

30、for(i=1;in;i+)vi=xi;if(norm=0)bet=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.42.2.4 例题例题把向量正交化为与平等且等长的向量的 Householder 变换为：，且其中：（1，-0.824343，-0.206086，-0.7213，-0.412171），=0.829128，11.7047.2.32.3 矩阵

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

32、eholder 变换使得.令，则有对于一般的矩阵householder 变换，假定我们已经进行了，使得步，得到了其中是上三角矩阵。假定第步是：先用算法 32.1 确定 Householder 变换使得其中；然后，再计算。令则其中将是上三角矩阵。这样，从约化为上三角阵。现在记出发，依次进行次，我们就可则注意，这样得到的上三角矩阵的对角元均是非负的。下面考虑计算的 QR 分解的存储问题。当分解完成后，一般来说，和。通常并不是将，而对每个算出，而是只存放和就不再需要，便可用来存放构成它的个 Householder 矩阵有如下形式，我们只需要保存即可。注意到我们正好可以将于存储在的对角元以下位置上。例

33、如，对的问题，其存储方式如下：综合上面的讨论，可得如下算法。算法算法 3.3.13.3.1(计算 QR 分解:Householder 变换法)容易算出,这一算法的运算量为2.3.32.3.3算法程序算法程序2.3.3.12.3.3.1 参数说明参数说明A：矩阵，作为输入是待分解的矩阵，作为输出上三角部分存储，下三角部分顺序存储各 Householder 向量Bet：顺序存储各 Householder 变换的M 和 n：指定矩阵的大小2.3.3.22.3.3.2 程序程序void HouseholderQR(double*a,double*bet,int m,int n)/m*n 阶矩阵的 QR

34、 分解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=Householder(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.42.3.4 例题例题矩阵：1 1-22

35、 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.927250.607625-3.23431 1.25684各 House 变换的 bet 值0.622036 0.162714 0.77532.42.4 求超定线性方程组的正交化算法求超定线性方程组的正交化算法2.4.12.4.1功能功能该算法用来求解超定方程组的最小二乘解2.4.22.4.2算法概述算法概述设分块为并且令有线性无关的列，并且假定。现将那么由此即知，

36、是最小二乘问题(3.1.3)的解当且仅当是的解。这求样一来，最小二乘问题(3.1.3)的解就可很容易从上三角方程组得。综合上面的讨论，可得正交化方法的基本步骤为：计算的 QR 分解（3.3.2）；计算；求解上三角方程组.2.4.32.4.3算法程序算法程序2.4.3.12.4.3.1参数说明参数说明A：超定方程组系数矩阵 b：超定方程组右端向量 bet：Householder 变换顺序值2.4.3.22.4.3.2程序程序void LeastSquares(double*a,double*b,double*bet,int m,int n)/最小二乘正交化算法HouseholderQR(a,be

37、t,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.42.4.4 例题例题用正交化法求方程组的最小二乘解。解结程序处理后得到：增广矩阵正交化后的结果：(左下三角部分各 Gouse 变换的 v 向量)11.7047 2.13589 6.15138-0.824343 3.79973-5.56321-0.206086-0.148331 0.133683-0.7213 1.36781 2.48331-0.4

38、12171 0.541989 0.16311各 House 变换的 bet 值0.829128 0.627619最小二乘解：0.79272-1.46411第第 3 3 章章 线性方程组的迭代解法线性方程组的迭代解法本章收集的算法是常见的单步线性定常单步线性定常迭代法即形如其中称为迭代矩阵迭代矩阵，是常数项常数项，是初始向量初始向量迭代法就是给出了一种构造向量序列的方法，对算法的最基本的要求就是算法的相容性，即迭代方程组（解的在具体选用算法还需要进一步考虑算法的收敛性，即由算法构造的向量序列使是有极限的，而且极限向量恰是方程组的解向量 即存在）与原方程组（）是同，且.单步线性定常迭代法的收敛性充

39、分必要条件是线性方程组的古典迭代算法，通常基于系数矩阵的如下分裂：其中，3.13.1 雅可比雅可比(Jacobi)(Jacobi)迭代法迭代法3.1.13.1.1 功能功能用雅可比迭代法求解线性方程组的解 但算法并不是万能的，其收敛性是有条件的详见教材3.1.23.1.2 算法概述算法概述雅可比迭代方程为，或迭代矩阵和常数项如下：，雅可比迭代公式，分量形式显然，雅可比迭代要求线性方程组系数矩阵的对角线元素非零并且迭代过程中，与各分量的计算次序无关3.1.33.1.3 程序设计程序设计3.1.3.13.1.3.1 参数说明参数说明A：方程组系数矩阵；b：方程组右端常向量；x：迭代解向量；n：方程

40、个数；Max：最大迭代次数；eps：最大允许误差3.1.3.23.1.3.2 程序程序double Jacobi(double*A,double*b,double*x,int n,doubleMax,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=fabs(x

41、i-yi);for(i=0;in;i+)yi=xi;return count;3.1.43.1.4 例题例题求解线性方程组迭代次数:33近似解向量:x0=0.413793 x1=0.172414 x2=0.862069误差限:1e-0063.23.2 高斯赛德尔高斯赛德尔(Gauss-Seidel)(Gauss-Seidel)迭代法迭代法3.2.13.2.1 功能功能用高斯赛德尔迭代法求解线性方程组的解 但算法并不是万能的，其收敛性是有条件的详见教材3.2.23.2.2 算法概述算法概述高斯赛德尔迭代方程为，或高斯赛德尔迭代公式，分量形式显然，高斯赛德尔迭代同样要求线性方程组系数矩阵的对角线元

42、素非零 但高斯赛德尔迭代与分量的计算次序有关，只能从小到大顺序计算3.2.33.2.3 程序设计程序设计3.2.3.13.2.3.1 参数说明参数说明A：方程组系数矩阵；b：方程组右端常向量；x：迭代解向量；n：方程个数；Max：最大迭代次数；eps：最大允许误差3.2.3.23.2.3.2 程序程序double GaussSeidel(double*A,double*b,double*x,int n,doubleMax,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+

43、)y=xi;xi=bi;for(j=0;ji;j+)xi-=Aij*xj;for(j=i+1;jn;j+)xi-=Aij*xj;xi/=Aii;if(errfabs(xi-y)err=fabs(xi-y);return count;3.2.43.2.4 例题例题求解线性方程组迭代次数:17近似解向量:x0=0.413793 x1=0.172414x2=0.862069误差限:1e-0063.33.3 超松驰超松驰(Succesive Over-Relaxation)(Succesive Over-Relaxation)迭代法迭代法3.3.13.3.1 功能功能用超松驰迭代法求解线性方程组的解

44、但算法并不是万能的，其收敛性是有条件的详见教材3.3.23.3.2 算法概述算法概述超松驰迭代方程为，或超松驰迭代公式，分量形式显然，同高斯赛德尔迭代一样，超松驰迭代要求线性方程组系数矩阵的对角线元素非零，而且与分量的计算次序有关，只能从小到大顺序计算3.3.33.3.3 程序设计程序设计3.3.3.13.3.3.1 参数说明参数说明A：方程组系数矩阵；b：方程组右端常向量；x：迭代解向量；n：方程个数；Max：最大迭代次数；eps：最大允许误差3.3.3.23.3.3.2 程序程序double SOR(double*A,double*b,double*x,int n,double w,dou

45、bleMax,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+)y=xi;xi=bi;for(j=0;ji;j+)xi-=Aij*xj;for(j=i+1;jn;j+)xi-=Aij*xj;xi=w*xi/Aii;xi+=(1.0-w)*y;if(errfabs(xi-y)err=fabs(xi-y);return count;3.3.43.3.4 例题例题求解线性方程组松驰因子:1.2迭代次数:10近似解向量:x0=0.413793 x1=0.172414 x2=0.

46、862069误差限:1e-006第第 4 4 章章 实矩阵特征值问题计算实矩阵特征值问题计算4.14.1 实矩阵上实矩阵上 H H 化的初等相似变换法化的初等相似变换法4.1.14.1.1 功能功能用初等相似变换把实矩阵化为上 Hessenberg 形矩阵（简称上 H-阵）。4.1.24.1.2 算法概要算法概要该方法需要对阶矩阵进行讨论算法的第一步，剩余的要以下几步：-2 次正交变换。由递归思想，我们只需要-1 步完全与之雷同。实现一次正交变换需选主元：寻找下标使，构造置换变换为，其中的第一个分量是所有分量中绝对值（或模）最大的。构造 Gauss 变换阵:,，其中，变换为。也就是说，经一次正

47、交约化后，将矩阵的第一列变换成上 H-阵列。4.1.34.1.3 算法设计算法设计算法 1（初等相似变换化 A 为上 H-阵）4.1.44.1.4 算法的算法的 C C 程序程序函数原型：void ElemHessenberg(double*a,int n);参数说明：a 既是输入矩阵，也是输出矩阵；n 是矩阵的阶数功能说明：运用初等相似变换将输入矩阵正交约化为上阵结果仍然保存在 a 中void ElemHessenberg(double*a,int n)int i,j,k;int p;double max,*l=new doublen;for(k=0;kn-2;k+)max=fabs(ak+

48、1k);p=k+1;for(i=k+2;in;i+)if(maxfabs(aik)max=fabs(aik);p=i;if(p!=k+1)for(j=k;jn;j+)max=apj;apj=ak+1j;ak+1j=max;for(i=0;in;i+)max=aip;aip=aik+1;aik+1=max;for(i=k+2;in;i+)li=-aik/ak+1k;for(j=k;jn;j+)aij+=ak+1j*li;for(j=k+2;jn;j+)for(i=0;in;i+)aik+1-=lj*aij;4.1.54.1.5 例题例题运用“初等相似变换法”化为上 H-阵：4.24.2 实矩阵上

49、实矩阵上 H H 化的列主序化的列主序 HouseholderHouseholder变换法变换法4.2.14.2.1 功能功能算法运用 Householder 变换从左到右逐列将阶实矩阵化为上 H-阵。4.2.24.2.2 算法概述算法概述设，Step1:构造-1 阶 Householder 变换（是，使-1 阶单位矩阵的第一列），并令，则至此，我们已经完成了第一步变换。并且有，值得注意的是：，就是说正交变换，不改变的 1 阶主子矩阵。Step2:构造-2 阶 Householder 变换（是，使-2 阶单位矩阵的第一列），并令，则至此，已完成了第 2 次正交变换。并且得到这里一直继续下去，直

50、到完成第。-2 步正交变换后，我们必将得到：。4.2.34.2.3 算法设计算法设计结合 Householder 变换构造算法，我们得到：4.2.44.2.4 算法算法 C C 程序程序函数原型：Void BeforeHessenberg(double*a,int n);参数说明：a 既是输入矩阵，也是输出矩阵；n 是矩阵的阶数功能说明：运用列序法将输入矩阵正交约化为上阵结果仍然保存在动态数组 a 中void BeforeHessenberg(double*a,int n)double*x=new doublen;double*v=new doublen;double b,t;int i,j,