c语言编程求解线性方程组论文

1、数值计算小报告题目：线性方程组求解方法比较 姓 名 和 丽 专 业 软件工程 班 级 11级软件（2）班 完成日期：2013 年 5月 18日摘 要 目前在许多实际应用领域,诸如航空、造船以及其它结构工程中,常遇到求解大型线性代数方程组的问题。本文根据线性代数方程组的雅可比迭代法、LU分解法及高斯列主元消去法三种解法进行了比较，用以方便在实际生活应用中更好的作出选择。在第二章中本文详细的介绍了线性代数方程组的三种解法的理论知识与证明过程。为了更加清晰的展现三种方法的不同点以及其各自的优越性，本文在第三章中给出了实例，通过实例的计算与程序的实现，再结合三种方法的优缺点进行了比较。 关键字：线性代

2、数方程组、迭代法、LU分解法、高斯列主元消去法、不同点、比较 目 录第一章 绪论 .4第二章 求解线性方程组的基本理论2.1 迭代法 .5 2.2 直接三角分解法 .6 2.3 高斯消去法 .7 第三章 三种算法求解方程组实例3.1 迭代法 8 3.2 直接三角分解法 .10 3.3 高斯列主元消去法 .14 3.4 三种方法的优缺点比较 .16参考文献 .17计算机专业学年论文线性方程组求解方法比较第一章 绪 论计算数学是数学学科的一大分支，它研究如何借助于计算机求解各类数值问题。应用计算机求解各类数值问题需要经历以下几个主要过程：1、实际问题2、数学模型3、计算方法4、算法设计5、计算求解

3、目前已有的数学软件可以帮助我们实现上机计算，基本上已经将数值分析的主要内容设计成简单的函数，只要调用这些函数进行运算便可得到数值结果。数值分析的内通包括线性代数方程组求解、非线性代数方程（组）求解、矩阵的特征值与特征值向量的计算、函数插值、函数逼近、数值积分与数值微分以及微分方程数值解法。线性方程组的求解从理论上可分为两类：直接法和迭代法。直接法是不考虑计算过程中的舍入误差，经过有限次的运算得到方程组精确解的方法，常见的方法是高斯顺序消去法、高斯列主元消去法和矩阵的LU分解法。迭代法是采用某种极限过程，用线性代数方程组的近似解逐步逼近精确解的方法。迭代法中常见的方法有简单迭代法、J-迭代法、G

4、S-迭代法和SOR-迭代法。本文主要是分析高斯列主元消去法、矩阵的LU分解法和简单迭代法理论上的异同，并用C语言程序通过具体实例进行了分析比较。本文将线性方程组的求解过程用计算机实现，本文的编写由以下几个特点：1、对于难点问题从具体模型引入，淡化抽象的概念与定理，通俗易通；2、对于具体模型本文给出了多种解题的思想及方法；3、对问题进行简洁易懂的理论证明，突出了线性代数的理论和基本思想，使数学方法更加利于理解掌握。4、简要分析了算法的计算效果、稳定性、收敛效果、计算精度以及优劣性。第1章 求解线性方程组的基本理论1.1 迭代法迭代法的基本思想：是将线性方程组转化为便于迭代的等价方程组，对任选一组

5、初始值（i=1,2n），按某种计算规则，不断地对所得到的值进行修正，最终获得满足精度要求的方程组的近似解。对于线性方程组Ax=b 其中，A为非奇异矩阵。将A分裂为A=M-N，其中，M为非奇异矩阵，且要求线性代数方程组Mx=d容易求解，一般选择为A的某一部分元素构成的矩阵，称M为A的分裂矩阵。于是，求解Ax=b转化为求解Mx=Nx+b,由此可构造一个迭代法： x(0)(初始向量) ， x(k+1)=Bx(k)+f (k=0,1,2) 其中，f=b/M,B=I-A/M为迭代法的迭代矩阵。选取M为A的对角元素组成的矩阵，即选取M=D，可得到解Ax=b的迭代法：x(0)(初始向量),x(k+1)=Bx

6、(k)+f (k=0,1,2) BJ为求解Ax=b的雅克比迭代法的迭代矩阵。解迭代法的计算公式为： （k=0,1,2，:i=1,2,3,.n）迭代法方法是求对称矩阵的全部特征值以及相应的特征向量的一种方法，它是基于以下两个结论：1)任何实对称矩阵A可以通过正交相似变换成对角型，即存在正交矩阵Q，使得 AQ=diag() 其中i（i=1,2,n）是A的特征值，Q中各列为相应的特征向量。2）在正交相似变换下，矩阵元素的平方和不变。即设，Q为交矩阵，记B=AQ=,则基本思想：是通过一次正交变换，将A中的一对非0的非对角线化成0，并且使得非对角元素的平方和减小。反复进行上述过程，使变换后的矩阵的非对角

7、元素的平方和趋于0，从而使该矩阵近似为对角矩阵，得到全部特征值和特征向量。 1.2 直接三角分解法 矩阵直接三角分解法：是高斯消去法的变形方法。高斯消去法有多种变形，有的是高斯消去法的改进，有的是用于某种特殊系数矩阵的化简。高斯消去法解线性方程组先消元，然后再回代。当用矩阵描述时，是对系数矩阵分解为一个上三角阵和一个下三角阵的乘积，即LU分解。因此，高斯消去法与矩阵的LU分解是一致的。将高斯消去法改写为紧凑形式，可以直接从矩阵A的元素得到计算L，U元素的递推公式，而不需要任何中间步骤，这就是所谓的直接三角分解法，一旦实现了矩阵A的LU分解，那么求解Ax=b的问题就等价于求解两个三角形方程组：

8、的问题，而这两个线性代数方程组只要回代，就可以求出其解。设A为非奇异矩阵，且有分解式A=LU，其中L为单位下三角矩阵，U为单位上三角矩阵，A= 第一步，用L的第一行分别乘以U的第j(j=1,2,n)列，比较两边可得(j=1,2,n) 分别用L的第i(i=1,2,n)行乘U的第一列，比较可得(i=1,2,n) 即得 (i=1,2,n) 这样就求出了L的第一列和U的第一行的所有元素。依次进行下去，一直到第k-1步，即已求出L的前k-1列和U的前k-1行的所有元素。第k步，用L的第k行分别乘U的第j(j=k,k+1,n)列，比较两边可得：=+ (j=k,k+1,n) 分别用L的第i(i=k+1,n)

9、行乘U的第k列，比较两边可得：=+ (i=k+1,n)总结上述讨论，得到用直接三角分解法求解Ax=b的计算公式： (j=1,2,n), (i=2,3,n) (j=k,k+1,n) (i=k+1,k+2,n)及求解Ly=b,Ux=y的计算公式： (k=2,3,n)， (k=n-1,n-2,2,1)1.3 高斯列主元消去法 高斯顺序消去法的基本思想是：对线性代数方程组所对应的增广矩阵（A|b）进行一系列“把某一行的非零常数倍加到另一行上”的初等变换，使得（A|b）中A的对角线一下的元素全变为0，从而使原方程组等价的转化为容易求解的上三角形线性代数方程组，再通过回代得到上三角形线性代数方程组的解，即

10、可求得原方程组的解。设线性方程组的增广矩阵为=首先，在第一列中选取绝对值最大的元素作为第一列的主元，即 然后交换第一行与第i行，经一次消元计算得：=(AB) ，此消去过程与高斯顺序消去法完全相同。 重复上述过程，设已完成第k-1步的选主元素，交换两行及消元过程后(AB)已约化为 第k步选主元素，在右下角方阵的第一列内选取绝对值最大的元素作为这一列的主元,即 = 然后交换的第i行与第k行，再进行消元计算。如此重复，直到最后将原线性代数方程组化为 = 回代求解得到 (k=n-1,2,1)列主元消去法除了每步需要按列选出主元，然后进行对换外，其消去过程与高斯顺序消去法是相同的。第2章 三种方法求解方

11、程组实例线性代数方程组的应用十分广泛，在实际生活中遇到的一些问题通常可以用求解方程组的方法来解决。本章主要是通过实例介绍三种方法的应用以及其各种方法的优缺点比较。例：用三种方法求解方程组 的解。2.1 迭代法 解题步骤： 将方程组记为Ax=b，其中 A= b= 将原方程组改写为 (1)也可写为 x=Bx+f 其中 B= f=任取初始值=,代入（1）式右边，得到新的值： 再将代入（1）式右边得到。反复利用这个计算程序，得到一个向量序列和一般的计算公式： =B+f 其中k表示迭代次数（k=1,2,）。 令，由雅克比迭代公式(k=0,1,2,)有： D 于是，根据雅克比迭代法的计算公式可求出方程组的

12、解。迭代法代码：#include<math.h>#include<stdio.h>int gausdl(int n,double i,double j,double x,double eps) int u,v,w,r; double p,t,s,q; for(u=0;u<=n-1;u+) w=u*n+u; p=0.0; xu=0.0;      for (v=0; v<=n-1; v+)      &

13、#160; if(u!=v)         r=u*n+v; p=p+fabs(ir);      if(p>=fabs(iw)        printf("failn"); return(-1);         p=eps+1.0;   

14、0;while (p>=eps)     p=0.0;      for (u=0; u<=n-1; u+)       t=xu; s=0.0;        for (v=0;v<=n-1;v+)       

15、;  if (v!=u) s=s+iu*n+v*xv;        xu=(ju-s)/iu*n+u;        q=fabs(xu-t)/(1.0+fabs(xu);        if (q>p) p=q;   return(1); void 

16、main() int u; double eps; static double  i16=2.0,0.0,0.0,1.0,2.0,3.0,2.0,4.0,-3.0,0.0,1.0, 6.0,1.0,-6.0,-5.0  static double x5,j4=0.0,-2.0,-7.0,6.0; eps=0.000001;if (gausdl(4,i,j,x,eps)>0)  for (u=0;u<=3;u+)   prin

17、tf("x(%d)=%13.7en",u,xu);运行结果： x1=-3.500000e+000x2=0.000000e+000x3=-3.000000e+000x4=7.000000e+000 2.2 直接三角分解法 解题步骤： 将方程组改写为= 设= 第一步有： 第二步有： 第三步有： 第四步有： 于是有 L= U= 由Ly=b得y，由Ux=y得x=LU分解法代码：Package com.nc4nr.chapter02.lu;public class LU .double a = .2.0, 0.0, 0.0, 1.0, 2.0, 2.0, 3.0, 2.0, 6.0,

18、 1.0, -6.0, -5.0, 4.0, -3.0, 1.0, 1.0 ;double b = 0.0,-2.0,6.0，-7.0； int anrow = 4; int indx = new intanrow; int parity = 1;private void lucmp() final double tiny = 1.0e-20; int imax = 0, n = anrow; double big, dum, sum, temp; double vv = new doublen; System.out.println("Origin coefficient matr

19、ix:"); output(a,4); for (int i = 0; i < n; i+) big = 0.0; for (int j = 0; j < n; j+) if (temp = Math.abs(aij) > big) big = temp; if (big = 0.0) System.out.println("lu: singular matrix in lucmp."); vvi = 1.0/big; for (int j = 0; j < n; j+) for (int i = 0; i < j; i+) sum

20、= aij; for (int k = 0; k < i; k+) sum -= aik*akj; aij = sum; big = 0.0; for (int i = j; i < n; i+) sum = aij; for (int k = 0; k < j; k+) sum -= aik*akj; aij = sum; if (dum = vvi*Math.abs(sum) >= big) big = dum; imax = i; if (j != imax) for(int k = 0; k < n; k+) dum = aimaxk; aimaxk =

21、ajk; ajk = dum; parity = -parity; dum = vvimax; vvimax = vvj; vvj = dum; indxj = imax; if (ajj = 0.0) ajj = tiny; if (j != n - 1) dum = 1.0/ajj; for (int i = j+1; i < n; i+) aij *= dum; output(a,4); private void lubksb() double sum; int n = anrow, ii = 0; System.out.println("Origin left-hand

22、 vector b:"); output(b,4); for (int i = 0; i < n; i+) int ip = indxi; sum = bip; bip = bi; if (ii != 0) for (int j = ii - 1; j < i; j+) sum -= aij*bj; else if (sum != 0.0) ii = i + 1; bi = sum; for (int i = n-1; i >= 0; i-) sum = bi; for(int j = i + 1; j < n; j+) sum -= aij*bj; bi =

23、 sum / aii; System.out.println("Final solution vector:"); output(b,4); private void output(double a, int anrow) for (int i = 0; i < anrow; i+) System.out.println(" | " + ai0 + " " + ai1 + " " + ai2 + " " + ai3 + " | "); System.out.printl

24、n("-"); private void output(double b, int bnrow) for (int i = 0; i < bnrow; i+) System.out.println(" | " + bi + " | "); System.out.println("-"); public LU() lucmp(); lubksb(); public static void main(String args) new LU(); 运行结果：Origin coefficient matrix: |

25、2.0 0.0 0.0 1.0 | | 2.0 2.0 3.0 2.0 | | 6.0 1.0 -6.0 -5.0 | | 4.0 -3.0 0.0 1.0 | -Origin left-hand vector b: | 0.0 | | -2.0 | | 6.0 | | -7.0 | -Final solution vector: | -3.5000000000000003 | | 0.0000000000000002 | | -3.0000000000000007 | | 7.0000000000000004 | - 2.3 高斯列主元消去法解题步骤：对线性方程组： Ax=b ； 令增广矩阵

26、为：=()消元过程： 令= ，把式的增广矩阵中的第1行的倍依次加到该增广矩阵的第i（i>1）行，则第i(i>1)行第j列位置的元素为： （i，j=2,3,4,n） (i=2,3,n)因此，式可转化为： = 依次做下去，一直做到第n-1步，即有 化为回代求解： （k=n-1，n-2，.,2,1）高斯列主元消去法代码：#include "stdlib.h"  #include "math.h"  #include "stdio.h"  

27、 int rgauss( int n,double a,double b)    int *js,l,k,i,j,is,p,q;     double d,t;     js=(int*)malloc(n*sizeof(int);     l=1;     for (k=0;k<=n-2;k+) 

28、60;    d=0.0;        for(i=k;i<=n-1;i+)        for(j=k;j<=n-1;j+)         t=fabs(ai*n+j);          if(t>d)d=t;jsk=j;

29、is=i;           if(d+1.0=1.0)l=0;    else      if(jsk!=k)       for(i=0;i<=n-1;i+)         p=i*n+k;q=i*n+jsk;    

30、;     t=ap;ap=aq;aq=t;             if(is!=k)       for(j=k;j<=n-1;j+)        p=k*n+j;q=is*n+j;        t=ap;ap=aq

31、;aq=t;       t=bk; bk=bis; bis=t;                if(l=0)      free(js);printf("failn");        retur

32、n(0);           d=ak*n+k;     for(j=k+1;j<=n-1;j+)       p=k*n+j;ap=ap/d;     bk=bk/d;     for(i=k+1;i<=n-1;i+)     

33、60; for(j=k+1;j<=n-1;j+)         p=i*n+j;         ap=ap-ai*n+k*ak*n+j;               bi=bi-ai*n+k*bk;       

34、0;   d=a(n-1)*n+n-1;    if(fabs(d)+1.0=1.0)     free(js);printf("failn");       return(0);         bn-1=bn-1/d;    for(i=n-2;i>=0;i-) 

35、60;    t=0.0;       for(j=i+1;j<=n-1;j+)         t=t+ai*n+j*bj;       bi=bi-t;          jsn-1=n-1;    for(k=n-1;k>

36、=0;k-)      if(jsk!=k)        t=bk;bk=bjsk;bjsk=t;    free(js);    return(1);  void main()   int i;    static double a16= 2.0,0.0,0.0,1.0,2.0,2.0,3.0,2.0,4.0,-3.0,0.0,1.0,6.0,1.0,-6.0,-5.0     static double b4=0.0,-2.0,-7.0,6.0 ;