数值线性代数课程设计_代数

《数值代数》课程设计评分标准 (1)交作业的内容 (1)附录1 论文结构撰写规范 (2)附录2 (2)参考论文1 (2)参考论文2 (13)参考论文3 (16)参考论文4 (21)1. 2-3两天查资料;2. 1-2天论文构思,列出提纲;3. 2-3确定计算方法和计算程序的详细内容和设计步骤；写课程设计报告；论文内容: 问题背景介绍, 建立正确的数学模型, 求解模型的数学原理, 计算过程,编写计算程序, 解释计算结果.论文字数: 1800-2200字. 注: 字数是指正文的字数.标题、摘要、关键字、目录、图表说明、参考文献、程序等不算入字数统计.论文格式: 5号字体,A4版面, 上下左右各留2cm,单倍行距,页面设置选为无网格; 编页码;图形尺寸：长度不超过4cm.一行可以放若干张图. 文稿中数据保留3位小数点.（1）要有勤于思考、刻苦钻研的学习精神和严肃认真、一丝不苟、有错必改、精益求精的工作态度，对有抄袭他人设计论文或找他人代设计、代做论文等行为的弄虚作假者一律按不及格记成绩。

（2）要敢于创新，勇于实践，注意培养创新意识。

（3）掌握课程的基本理论和基本知识，概念清楚，设计计算正确，结构设计合理，实验数据可靠，软件程序运行良好，绘图符合标准，论文撰写规范，答辩中回答问题正确。

附录1 论文结构撰写规范（1）题目、院系、班级、学生姓名。

（2）摘要摘要是论文内容的简短陈述，一般不超过150字。

关键词应为反映论文主题内容的通用技术词汇，一般为4个左右，一定要在摘要中出现。

（3）正文正文应按目录中编排的章节依次撰写，要求计算正确，论述清楚，文字简练通顺，插图简明，书写整洁。

文中图、表及公式不能徒手绘制和书写。

（4）参考文献（资料）参考文献必须是学生在课程设计中真正阅读过和运用过的，文献按照在正文中的出现顺序排列。

各类文献的书写格式如下：a．图书类的参考文献序号作者名·书名·（版次）·出版单位，出版年：引用部分起止页码。

b．翻译图书类的参考文献序号作者名·书名·译者·（版次）出版单位，出版年：引用部分起止页码。

c．期刊类的参考文献序号作者名·文集名·期刊名·年，卷（期）：引用部分起止页码。

附录2数值代数课程设计题参考论文1《数值线性代数》课程设计信息与计算科学0440501106 蔡云摘要：Gauss消去法是目前求解中小规模线性方程组（即阶数不要太高，例如不超过1000）常用的方法，它一般用于系数矩阵稠密（即矩阵的绝大多数元素都是非零的）而没有任何特殊结构的线性方程组。

选主元Gauss消去法弥补了不选主元的Gaauss消去法的不足，但是也付出了昂贵的代价。

平方根是求解对称正定线性方程组最常用的方法之一。

关键词：Gauss 消去法，选主元Gauss 消去法，平方根法正文：一，问题背景：在古代数学和近代数学数值计算和工程应用中，求解线性方程组都是重要的课题。

西方的Gauss 消元法在求解未知数多的大型线形方程组则是在20世纪中叶电子计算机问世后才成为现实。

这次实验将充分运用各种数值计算方法，及计算机的强大数值计算功能来解决所遇到的问题。

二，问题提出：1.先用你所熟悉的计算机语言讲不选主元和列主元Gauss消去法编写成通用的子程序，然后用你所编写的程序求解下面的84阶方程组最后，将你的计算结果与方程组的精确解进行比较，并就次谈谈你对Gauss 消去法的看法。

2．先用你所熟悉的计算机语言将平方根法和改进的平方根法编写成通用的子程序，然后用你编写的程序求解对称正定方程组Ax=b,其中 （1）b 随机的选取，系数矩阵为100阶矩阵1011101110111011101110⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦；（2）系数矩阵为40阶Hilbert 矩阵，即系数矩阵A 的第i 行第j 列元素为11i j a i j =+-，向量的第i 个分量为111ni j b i j ==+-∑。

3．用第1题的程序求解第2题的两个方程组并比较所有的计算结果，然后评论各个方法的优劣。

三，模型建立与数学原理： （1）Gauss 消去法Gauss 消去法即是把一个给定的矩阵分解为一个下三角L 和一个上三角U 的乘积。

为此我们找到了一个最简单的下三角矩阵T k k k L I l e =-.其中1,(0,,0,,,)k k k nk l l l =+即1,,1111k k kn k L l l +⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥-⎢⎥⎢⎥⎢⎥-⎢⎥⎣⎦这种类型的下三角称作Gauss 变换。

而称向量k l 为Gauss 向量。

对于一个给定的向量1(,,)T n n x x x R =∈ ，我们有111,(,,,,,)T k k k k k k n k nk L x x x x x l x x l ++=-- 只要取,1,,i ik k xl i k n x ==+ ，便有1(,,,0,,0)T k k L x x x = 。

于是将Gauss 变换作用于一个矩阵，就将对应的列向量相应项变为0。

一般的阶矩阵A,在一定条件下，我们也可以计算n-1个Gauss 变123828384617861158611586115861158614x x x x x x ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦换121n L L L - 使121n n L L L A -- 为上三角矩阵。

A LU = ，这种计算主元分解的方法称做Gauss 消去法。

LU 分解的可行性分析：若n n A R ⨯∈的顺序主子阵L 均非奇异，则存在唯一的单位下三角n nL R⨯∈和上三角阵n nU R⨯∈使得A LU =对于题1中，我们可以通过高等代数的方法计算出对于*det()A =*.det()A ≠0,也即存在唯一性。

(2) 列主元Gauss 消去法由于可能存在小主元现象，在Gauss 消去的过程除数的绝对值过小，则导致其他变量量级的巨大增长和舍入误差扩散。

因而要保持数值解得精度应避免选取绝对值过小的元素作主元。

(3) 平方根法Cholesky 分解定理：若n nA R ⨯∈对称正定，则存在一个对角元均为正数的下三角阵n nL R ⨯∈，使得TA LL =。

Cholesky 分解可用不选主元的Gauss 消去法来实现，更简单的方法是比较A=LL T 两边对应元素来计算L 的，设11212212l n n nn l l l L l l ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦比较TA L L =两边对应元素,得关系式1nij ip jpp a ll ==∑ 1j i n ≤≤≤。

经计算得11()/k ik ik ip kp kk p l a l l l -==-∑ ，1,,i k n =+这样便求出了L 的第k 列元素。

这种方法称为平方根法。

四，计算过程及敏度分析： (1)在1中对题目用高斯消去法C++编程（见附录一）后运行，结果用matlab 作图后显示可以看出，在初始分量，效果还是比较理想的逼近于x=1,到后面出现了震荡。

可见得出的解数值不太稳定。

这是由于系数矩阵的条件数κ(A)非常大，用matlab 计算（程序见附录六）可得κ(A)=cond(A)= 4.1779e+016，已经非常大，也即该线形方程组得求解问题是病态的，我们得到的数值解是不稳定的。

而应用列主元高斯消去法，当矩阵是84阶时，不能得出结果，这是因为条件数很大时矩阵接近于一个奇异矩阵，高斯消去法不一定能够进行下去。

当为30阶时，可以看出也得到了精确解（附录中只提供了84阶矩阵，相应30阶读者可以自行改动）。

(2)2中系数矩阵是严格对角占优的。

条件数κ(A)=Cond(A)=1.4997,比较小，可知此问题数值比较稳定。

由于细数矩阵严格对角的方程组Gauss-seidel 迭代收敛，所以可以应用此算法来解决（具体程序见附录）。

从计算结果的图形表示我们可以看出，此种方法达到了数值要求，比较好的收敛于线性方程的解。

应用平方根法时的程序见附录四。

由于b 是计算机随机选取的，所以得到的解是人为无法确定的，我们在这里只画出了一种解的情况下的图象（见图3）用改进的平方根法解Hilbert 矩阵的程序见附录五，相应的图形件图4。

有趣的是，当我们用平方根法解十并不会得到结果。

应着重指出，正定矩阵的TLL分解具有重要的性质，即ik l ≤L 的元素的绝对值不会变得很大，这个性质对计算方法的稳定性是非常重要的。

（3）平方根法可以免去高斯消元过程。

但是只限于对称方程组。

高斯消元法是比较经典的求解稠密线性方程组的方法。

为避免小主元，应选用高斯列主元法。

五，参考文献：[1] 徐树方等·数值线性代数·北京大学出版社·2000年;[2]张韵华，陈效群等·数值计算方法和算法·科学出版社·2002年[3] 蒋长锦·科学计算和C程序集·中国科学技术大学出版社·1998年六，附录：附录一：Gauss消去法#include<iostream.h>#include<stdlib.h>#include<math.h>#ifndef Epsilon#define Epsilon 0.00000001#endif#define DIM 84int gcplim(int process,double A[DIM][DIM],double xx[DIM]){int k,i,j;//double pelement;if(process==1)cout<<"The process of elimination:\n";for(k=0;k<DIM;k++){/*pelement=fabs(A[k][k]);i0=k;for(i=k;i<DIM;i++)if(fabs(A[i][k])>pelement){pelement=fabs(A[i][k]);i0=i;}if(i0!=k){for(j=0;j<DIM;j++){pelement=A[k][j];A[k][j]=A[i0][j];A[i0][j]=pelement;}pelement=xx[k];xx[k]=xx[i0];xx[i0]=pelement;}*/if(fabs(A[k][k])<Epsilon) return(1);for(i=k+1;i<DIM;i++){A[i][k]=A[i][k]/A[k][k];for(j=k+1;j<DIM;j++)A[i][j]=A[i][j]-A[i][k]*A[k][j];xx[i]=xx[i]-A[i][k]*xx[k];}if(process==1){for(i=0;i<DIM;i++){for(j=0;j<DIM;j++)cout<<A[i][j]<<'\t';cout<<" |"<<xx[i];cout<<endl;}cout<<endl;}}//for(i=DIM-1;i>=0;i--)//{xx[DIM-1]=xx[DIM-1]/A[DIM-1][DIM-1];for(i=DIM-2;i>=0;i--){for(j=i+1;j<DIM;j++)xx[i]=xx[i]-A[i][j]*xx[j];xx[i]=xx[i]/A[i][i];}return(0);}void main(){int i;static double A[DIM][DIM];for(i=1;i<DIM;i++){A[i][i]=6;A[i][i+1]=1;A[i][i-1]=8;}A[0][0]=6;A[0][1]=1;A[DIM-1][DIM-1]=6;A[DIM-1][DIM-2]=8;static double b[DIM];for(i=1;i<DIM-1;i++)b[i]=15;b[0]=7;b[DIM-1]=14;cout<<"Colume"<<endl;if(gcplim(1,A,b)==1){cout<<"The linear system hasn't solution!"<<endl;cout<<"Strike any key to exit!"<<endl;exit(1);}cout<<"Column principle elimination for the solution:\n";for(i=0;i<DIM;i++)cout<<b[i]<<endl;}解的图象如下所示：图表1附录二：Gauss列主元消去法#include<iostream.h>#include<stdlib.h>#include<math.h>#ifndef Epsilon#define Epsilon 0.00000001#endif#define DIM 84int gcplim(int process,double A[DIM][DIM],double xx[DIM]){int k,i,j,i0;double pelement;if(process==1)cout<<"The process of elimination:\n";for(k=0;k<DIM;k++){pelement=fabs(A[k][k]);i0=k;for(i=k;i<DIM;i++)if(fabs(A[i][k])>pelement){pelement=fabs(A[i][k]);i0=i;}if(i0!=k){for(j=0;j<DIM;j++){pelement=A[k][j];A[k][j]=A[i0][j];A[i0][j]=pelement;}pelement=xx[k];xx[k]=xx[i0];xx[i0]=pelement;}if(fabs(A[k][k])<Epsilon) return(1);for(i=k+1;i<DIM;i++){A[i][k]=A[i][k]/A[k][k];for(j=k+1;j<DIM;j++)A[i][j]=A[i][j]-A[i][k]*A[k][j];xx[i]=xx[i]-A[i][k]*xx[k];}if(process==1){for(i=0;i<DIM;i++){for(j=0;j<DIM;j++)cout<<A[i][j]<<'\t';cout<<" |"<<xx[i];cout<<endl;}cout<<endl;}}for(i=DIM-1;i>=0;i--){//xx[DIM-1]=xx[DIM-1]/A[DIM-1][DIM-1];//for(i=DIM-2;i>=0;i--)// {for(j=i+1;j<DIM;j++)xx[i]=xx[i]-A[i][j]*xx[j];xx[i]=xx[i]/A[i][i];}return(0);}void main(){int i;static double A[DIM][DIM];for(i=1;i<DIM;i++){A[i][i]=6;A[i][i+1]=1;A[i][i-1]=8;}A[0][0]=6;A[0][1]=1;A[DIM-1][DIM-1]=6;A[DIM-1][DIM-2]=8;static double b[DIM];for(i=1;i<DIM-1;i++)b[i]=15;b[0]=7;b[DIM-1]=14;cout<<"Colume"<<endl;if(gcplim(1,A,b)==1){cout<<"The linear system hasn't solution!"<<endl;cout<<"Strike any key to exit!"<<endl;exit(1);}cout<<"Column principle elimination for the solution:\n";for(i=0;i<DIM;i++)cout<<b[i]<<endl;}当DIM=84时程序最后输出结果为：The linear system hasn't solution!Strike any key to exit!Press any key to continue附录三：Gauss-seidel迭代#include<iostream.h>#include<math.h>#include<stdlib.h>#define DIM 100#define E 1E-6double f(double x[DIM],double y[DIM]){double m[DIM];for(int i=0;i<DIM;i++){m[i]=fabs(x[i]-y[i]);if(m[i]>m[0])m[0]=m[i];}return m[0];}int jacobi(double a[DIM][DIM],double b[DIM],double x[DIM]){static int first=0;int i,j,k;static double x1[DIM],norm;norm=0;if(first==0){for(i=0;i<DIM;i++){x[i]=rand();norm+=x[i]*x[i];}norm=sqrt(norm);for(i=0;i<DIM;i++)x[i]/=norm;first=1;}for(k=0;k<10000;k++){for(i=0;i<DIM;i++){x1[i]=b[i];for(j=0;j<i;j++)x1[i]-=a[i][j]*x[j];//若为Gauss-Seidel迭代，则此行变为"x1[i]-=a[i][j]*x1[j];"其他地方不变。