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

1．矩阵的三角分解里列主元三角分解：，这里位下三角矩阵且各元素的模不超过1，阵。，，2．矩阵的正交三角分解和一个上三角矩阵的把实矩阵分解成单位下三角形矩阵和上三角形矩阵．阶方阵作如下分解则角矩阵为止。即有．．．．综上所述，我们可综上所述，我们可以将两个矩阵因子继续存储在原矩阵的存储空1.1.3.1.1参数说明1.1.3.1.2C程序{rintkknk{falsekini{kkkjnjaijaikakj;}}ntrue}，这里是单位下三角矩阵且各元素的模不超阵。所不同的是在构造Gauss变换前，先在对应列中选择绝对值最大的元素(称为列主元)，然选主元：确定使；行交换：将矩阵的第行和第行上的元素互换位置。即．1.2.3.1参数说明1.2.3.2C程序{nipiiiini{maxaii];{maxaji];}urni{}iiajkaikaji;}returnn-1;rnn}P3)，，法基本上类似，所不同的是在构造Gauss变换前，先在对应列中选择绝对值最大的元素(称为列主元)，然使行列交换：将矩阵的第行和第行上的元素互换位置，再将第列和1.3.3.1参数说明1.3.3.2C程序intAllGaussLU(double**a,int*p,int*q,intn)//n阶矩阵的全主元{qiiiini{maxaii];iqii{maxajk];jqik}urni}}iiajkaikaji;}returnn-1;rnn}256811166.750.754440.7916670.7222220.6315790.2105261.4对称正定Cholesky分解，，，的元素(的元素(另一方面，则的元素(也可由Gauss消去法的思想逐列得到，1.4.3.1参数说明1.4.3.2C程序{for(k=0;k<n;k++){sqrtakkkkaijaikajk;}}．，设，以列为序从左到右，用待定系数法，则的元素(以及对角矩阵的元素可由下列公式确定：1.5.3.1参数说明1.5.3.2C程序{blennj{ijiviajiaiijjajivi{kijaikvkjj}}deletev}，，第第2步：再消去后个方程中的，即，，第，其增广矩阵为fork=1:n-1bibiai,k)b(k)voidLowerdoubleadoubleb,intn){orintiinibiaii;ijnjbjbiaj][i];}b[n-1]/=a[n-1][n-1];}voidUnitLowdoubleadoubleb,intn){iniijnjbjbiaj][i];}个方程中的，即第第，第第2步：再消去前个方程中的，即，个方程中的，即，，，其增广矩阵为角矩阵时，算法或简化为voidUperdoubleadoubleb,intn){forintin;i>0;i--){biaii;ijbjbiaj][i];}ba[0];}的回代法voidUnitUpdoubleadoubleb,intn){forintin;i>=0;i--)ijbjbiaj][i];}等于零。1.8.2.1算法基于系数矩阵的三角分解。将作三角分解：；；．1.8.3.1参数说明1.8.3.2C程序boolGaussdouble*a,double*b,intn){nnitLowabnperabnntrue}false}上三角方程组上三角方程组的解(即原方程组的解)：(-1，1，-1，1)．((1) 其中矩阵和的形式如下：，而和的元素由下列关系式确定：． (3)进而解方程组(1)就等价于解方程组： (6)式为“追”的过程，(7)式为“赶”的过程，故称该法为追赶法。1.9.3.1参数说明.3.2C程序voidPursuitdoubleadoublebintn/追赶法{wa1][0];ba[0];kknk{akwwa[1][k+1]-a[0][k]*a[1][k];bkbk1]-a[0][k]*b[k])/w;}forknk>=0;k--)bkbka][k]*b[k+1];}1.10列主元高斯消去法用列主元三角分解法(算法5_2)对系数矩阵进行列主元三角分解：即对向量进行行调整：；求解单位下三角方程组：；1.10.3.1参数说明1.10.3.2C程序boolRowGaussdoubleadouble*b,intn)//列主元高斯消去法{falserintiini{bibpiit}nitLowabnperabnntrue}，其中，；的解为：=(2，188，-38.5714，2.18182)；的解，即原方程组的解：(-1，1，-1，1)。用列主元三角分解法(算法5_2)对系数矩阵进行列主元三角分解：即对向量进行行调整：；求解单位下三角方程组：；1.11.3.1参数说明1.11.3.2C程序boolAllGaussdoubleadouble*b,intn)//全主元高斯消去法{falserintiini{{bibpiit}}nitLowabnperabnforini>=0;i--){{bibqiit}}ntrue}256811166.750.754440.7916670.7222220.6315790.210526线性方程组。eran1.12.3C程序voidCholeskydoubleadoubleb,intn){skyLLanbniniijnjiperabn}线性方程组。eranntn{yLDLannitLowabnini{biaii;ijnji}nitUpabn}，，2.1.1功能通过解法方程组的解求得方程组算法概述定理是的最小二乘解当且仅当算法1(正则化法：求解的最小二乘解)voidLeastSquaresdoubleadoubleb,double*x,intm,intn){oublenini{nj{xi0;k{iakjxiakibk}}}yLLcnxnijnjpercxncideletec}1.4例题--00-027550081252290000-00022.2Householder变换2.2.1功能针对非零向量，构造Householder变换，使算法概述设HouseholderHouseholder变换满足，，证明令则■定理证明表明，对任意的定理证明表明，对任意的可按如下的步骤来构造确定：计算计算时，为避免溢出，可以考虑先将规范化为使，h在计算向在计算向量的第一个分量时，为了避免两相近数相减，而，，算法程序2.2.3.1参数说明2.2.3.2C程序ouseholderdoublexdoublevintn{fabsxiabsxi{AfxMessageBox("Fail,becausexiszero.\n");}norm0;mxixiv.0;bet=0.0;{hsqrtxxnormxalphvnormxalph;betv*v[0]/(norm+v[0]*v[0]);lph}t}2.4例题2.3.1功能阵．即求得正交阵．即求得正交矩阵使算法概述定定理(QR分解定理)设且当且当且非奇异时，上述的分解还是唯一的。■..令对于一般的矩阵，假定我们已经进行了步，得到了householder变换，使得其中是上三角矩阵。假定则其中其中是上三角矩阵。这样，从出发，依次进行次，我们就可将将约化为上三角阵。现在记则阵的对角元均是非负的。就不再需就不再需要，便可用来存放和。通常并不是将算出，而是只存放和即即可。注意到有如下形式我们正好可以将我们正好可以将存储在的对角元以下位置上。例如，对算法程序2.3.3.1参数说明2.3.3.2C程序{ublemublem{kbetkHouseholderxk,v+k,m-k);{i}}}3.4例题10-1211.13389252.64575-2.22045e-0166250.9114381.92725234311.25684140.77532.4.1功能算法概述设这计计算的QR分解(3.3.2)；求求解上三角方程组.算法程序2.4.3.1参数说明er换换顺序值2.4.3.2C程序{HouseholderQRabetmnnj{jimibjtbitaij];}perabn}4.4例题433.79973-5.56321-0.1483310.133683.367812.483310.5419890.16311191的算法是常见的单步线性定常迭代法．即形如算法的相容性，即迭代方程组(算法的相容性，即迭代方程组()与原方程组()是同使..单步线性定常迭代单步线性定常迭代法的收敛性充分必要条件是．线性方程组的古典迭代算线性方程组的古典迭代算法，通常基于系数矩阵的如下分裂：．．，3.1雅可比(Jacobi)迭代法1.1功能．常数项，，3.1.3.1参数说明3.1.3.2C程序JacobidoubleAdoublebdoublexintndoubleMax,doubleeps){blenrriwhileerr=eps&&count<Max){{xiAijyjxiAijyjii}}count}例题x.413793x[1]=0.172414x[2]=0.8620693.2高斯－赛德尔(Gauss-Seidel)迭代法2.1功能，，阵的对角线元素非3.2.3.1参数说明3.2.3.2C程序GaussSeideldoubleAdoublebdoublexintndoubleMax,doubleeps){yerrxi0;whileerr=eps&&count<Max){{xiAijxjxiAijxjii}}count}例题x0]=0.413793x[1]=0.172414x[2]=0.8620693.3超松驰(SuccesiveOver-Relaxation)3.1功能，或．，3.3.3.1参数说明3.3.3.2C程序SORdoubleAdoublebdoublexintndoublewdoubleMax,doubleeps){yerrxi0;whileerr=eps&&count<Max){{xiAijxjxiAijxjxiwxiAiiiwy}}count}例题松驰因子:1.2x.413793x[1]=0.172414x[2]=0.8620694.1.1功能用初等相似变换把实矩阵化为上Hessenberg形矩阵(简称上H-阵)。4.1.2算法概要选主元：寻找下标选主元：寻找下标使，构造置换变变换为，，，其其中的第一个分量是所有分量中绝对值(或模)最大的。，，变变换为4.1.3算法设计ElemHessenbergdoubleaintnemaxlnewdoublenfor(k=0;k<n-2;k++){maxfabsak][k]);maxfabsaik]);}maxapj];akj}maxaip];aik}}for(i=k+2;i<n;i++){forjkjnjaijakj]*l[i];}设aikljaij;}}}4.1.5例题4.2.1功能4.2.2算法概述，，，r(是-1阶单位矩阵的第一列)，并令，，则且有，，，(是-2阶单位矩阵的第一列)，并令，，则这里。4.2.3算法设计voidBeforeHessenbergdoubleaintn{blenblenkknk{x[i-k]=a[i][k-1];bBeforeHousexvnk);{iniijtvik}aik.0;{jnjijtvjk}}}4.2.5例题4.3.1功能4.3.2算法概述Setp记是-1阶单位阵的第-1列。．Step-2，便逐步得到上H-阵。4.3.3算法设计voidBehindHessenbergdoubleaintn{wdoublenforknk>1;k--){bBehindHouseakvk{ikii}{jkjj}jjkjakj;}}4.4Householder变换程序BeforeHousedoublexdoublev,intn);BeforeHousedoublexdoublevintn{bsxiniinimxi*x[i];rtxxmvx0]-a;vmxabvv0]/(m+v[0]*v[0]);}b}BehindHousedoublexdoublev,intn);BehindHousedoubleydoublevintn{bleninibsxformiini+)mxi*x[i];vn1;iiniasqrtxnxn1]+m);vnx[n-1]-a;vnm/(x[n-1]+a);b=2.0*v[n-1]*v[n-1]/(m+v[n-1]*v[n-1]);v[i]/=v[n-1];}deletexb}4.5求实上Hessenberg矩阵特征值的QR4.5.1功能4.5.2算法概要 (1.1)与相似，亦即任一与原始矩阵若若原始矩阵是上H－矩阵，则容易证明：所有也是上H－矩阵．若矩阵是上H－阵，对它进行QR分解时，工作量是比较小的．因此，对于一般矩阵，应先把它化为上H－阵，然后再采用QR算法．或因为实矩阵因为实矩阵可能有复特征值，如果变换(1.1)带有复位移量(一般取取的一粗糙的近似特征值)，则一般是一复矩阵．在计算上的两步合成一步两步合成一步)中，使用了两个实的和或一对复共轭的值，而避两步QR算法以下述论据为基础．若矩阵B非奇异，且以确为上H－阵，则和次对角线元素，则确定是唯一的．现在叙述两步QR算法，基基本关系如下：,(1.5)，其中是正交阵，是上H－阵，于是，如果具有和相同的第一列(亦即的第一行与的第一行相同)，则根据上述论据，一定有＝＝及(亦即(亦即)．从(1.7)式可知，是实的或复共轭时，都是一个实矩阵．因且行即行即的第一行，因此，我们只需要导出满足方程(1.8)(或(1.9))且其第第一行元素与的第一行元素相同的矩阵．用以用以三角化任一矩阵的第一个因子仅由该矩阵(例如)的第一列所由于由于的第一列只有前三个元素非零，所以与相联系的也只有在这是一个H－阵再加三个元素(＂＋＂元素)的矩阵．又因为任一实矩阵又因为任一实矩阵可由一正交相似变换化为上H－阵，其形式如下：其中由(1.10)定义．特别，我们取，则有并且的第一行即的第一行．由此可知，就是满足方程(1.8)(或(1.9))的矩阵QR一定是．这就是关于Francis提出的两步QR方法．显然，该算法在任一阶段都不出现复数．为为上H－矩阵时与相联系的的确定方法．相似变换后，都产生具有下面形式的矩矩阵．例如以8阶矩阵为例，在进行变换前的矩阵形状为这是一个具有在(5，3)，(6，3)，(6，4)位置上的三个附加元素的据此，对于与Householder变换相联系的具有下面的形式：它只在上具有非零元素，这些元素可由中位于，关于原点位移量的选择，在每一阶段，是由关于原点位移量的选择，在每一阶段，是由右下角的2阶子矩阵的两可以被可以被忽略，则该迭代可认为已达目的．当可忽略，则可作为一特征值，并去掉最可作为一特征值，并去掉最后一行和最后一列而把矩阵收缩为阶上的的两个特征值．进而去掉最后两行或两列而把矩阵收缩为阶上H－看看是否可能被忽略(见说明)．若有，则矩阵的特征值问题可分裂为两素的乘积可能较小，以致允许我们在较低的子矩阵进行迭代(见说明)．在利用上面提到的在利用上面提到的QR算法时，和都将是零，因为该矩阵是正交的，用没有原点位移的QR算法，其结果将是一个不变式．式式来确定和.5.3算法主程序及使用说明4.5.3.1voidfrancis(double**a,double*wr,double*wi,int*cn,intn,doubleeps=1.0e-13)算得知的上算得知的上H－阵的特征值计算出来，并将第个特征值的实部和虚部分cn[i]中．4.5.3.2voidhqr(double**h,doubles,doublet,intn,intl,intm)4.5.3.3voidhessenberg(double**a,intn)算上H－阵的特征值，不需要调用该函数．4.5.3.4doublehouse(double*x,double*v,intn)．．4.5.3.5voideigenvalue2(double**h,double*wr,double*wi,intm)该函数负该函数负责计算出阶矩阵的两个特征值，并依次把它们实部放入4.5.3.6voidmain()4.5.4例题-1781278121078121-107812迭代次数03-32-23-177400000000迭代次数02-22-21-173535e-0060000迭代次数002-2，，75241541511311306955409106-65349434951305305691565235652532055305507139333936982323232127604323604323273221232215723372298538370.93366404470.39254638454845355767471999095389-43-34-45-53-33-3-18当当时01-1502950当当时迭迭001553-105时0739100155.5算法程序ostreamhomaniphstreamhathhvoidhqrdoublehdoublesdoubletintnintl,intm);ousedoublexdoublevintnvoidhessenbergdoubleaintnvoideigenvaluedoublehdoublewrdoublewi,intm);voidviewmdoubleaintn{oublenini{ublennj}berganublenublencnnforiinicoutsetprecisionsetw5)setwwiisetw}{doubles,t;//调用hqr()所需的参数，用以确定位移量mnl;whilem0){if(m==1)//此时待处理矩阵已压缩为一阶(即一个实数){wr[m-1]=h[m-1][m-1];wi[m-1]=0.0;m--;}elseifm矩阵，可直接计算特征值{mm=2;}elseifdhmmhmm]){forlml>0;l--)abshmmfabshmm}shmmhmm1];thmmhmm]-h[m-2][m-1]*h[m-1][m-2];}{}rhstnlmimibshiiepsfabshiifabshiihi[i-1]=0.0;}hmm{mm=2;}wr[m-1]=h[m-1][m-1];wi[m-1]=0.0;m--;}}}voidhqrdoublehdoublesdoubletintnintl,intm){ubleublexhllhl[l]+h[l][l+1]*h[l+1][l]-s*h[l][l]+t;x=h[l+1][l]*(h[l][l]+h[l+1][l+1]-s);xhllhl2][l+1];housexv{iliilhijh[i][j]-=a*v[i-l];}{jljahi[j]*v[j-l];h[i][j]-=a*v[j-l];}k{kixikhik];housexv{ikikiah[i][j]*v[i-k-1];kih[i][j]-=a*v[i-k-1];}{jkjkjah[i][j]*v[j-k-1];kjh[i][j]-=a*v[j-k-1];}hi[k]=0.0;}ximh[i][m-3];housexv{mimijvimhijavim2];}{mjmjjvjmhijavjm2];}hm1][m-3]=0.0;}housedoublexdoublevintn{bsxiniinimxi*x[i];rtxxmvx0]-a;vmxabvv0]/(m+v[0]*v[0]);}b}voidhessenbergdoubleaintn{blenblenkknk{x[i-k]=a[i][k-1];bhousexvnk);{iniijtvik}aik.0;{jnjijtvjk}}}voideigenvaluedoublehdoublewrdouble*wi,intm){shmmhmm1];thmmhmm]-h[m-2][m-1]*h[m-1][m-2];wi[m-1]=wi[m-2]=0.0;{wr[m-2]=(s+sqrt(d))/2.0;wr[m-1]=t/wr[m-2];}{wr[m-1]=(s-sqrt(d))/2.0;wr[m-2]=t/wr[m-1];}}wr[m-2]=wr[m-1]=s/2.0;wi[m-2]=sqrt(-d)/2.0;wi[m-1]=-wi[m-2];}}voidviewmdoubleaintn{ini{