高斯列主消元法法和LU分解法.doc
【精品文档】如有侵权,请联系网站删除,仅供学习与交流高斯列主消元法法和LU分解法.精品文档.高斯列主消元法:#include<stdio.h>#include<math.h>#define N 10int main() int i=0,j=0,r=0,n=0,k=1; float aNN,xN,max=-1,sum1=0,sum2=0;/输入矩阵维数 printf("Please enter n 0<n<11:n"); scanf("%d",&n);/输入矩阵元素,最后一列为b for(i=1;i<=n;i+) for(j=1;j<=n+1;j+) printf("Please enter a%d%d:n",i,j); scanf("%f",&aij); do for(i=k;i<=n;i+) if(fabs(aik)>max) max=aik; r=i; if(ark=0) printf("Cannot solve!n"); return 0;/换行 for(j=1;j<=n+1;j+) sum1=akj; akj=arj; arj=sum1;/将矩阵化为上三角矩阵 for(i=k+1;i<=n;i+) for(j=k+1;j<=n+1;j+) aij=aij-(aik)*(akj)/(akk); k=k+1; max=-1; while(k<(n-1);/求x xn=(ann+1)/(ann); for(k=n-1;k>=1;k-) for(j=k+1;j<=n;j+) sum2=sum2+(akj)*(xj); xk=(akn+1-sum2)/(akk); sum2=0;/输出结果 for(i=1;i<=n;i+) printf("x%d=%fn",i,xi); while(1); return 0;LU分解法:#include <stdio.h> #include <stdlib.h> #define N 10 /矩阵大小范围 * 使用已经求出的x,向前计算x(供getx()调用) * float a 矩阵U * float x 方程组解 * int i 解的序号(数组X元素序号) * int n 矩阵大小 * return 公式中需要的和 float getmx(float aNN, float xN, int i, int n) float mx = 0; int r; for(r=i+1; r<n; r+) mx += air * xr; return mx; * 使用已经求出的y,向前计算y(供gety()调用) * float a 矩阵L * float y 数组Y * int i 数组Y元素序号 * int n 矩阵大小 * return 公式中需要的和 float getmy(float aNN, float yN, int i, int n) float my = 0; int r; for(r=0; r<n; r+) if(i != r) my += air * yr; return my; * 解方程组,计算某x * float a 矩阵U * float x 方程组解 * int i 解的序号 * int n 矩阵大小 * return 方程组的第i个解(数组X的第i个元素值) float getx(float aNN, float bN, float xN, int i, int n) float result; if(i=n-1) /计算最后一个x的值 result = (float)(bi/an-1n-1); else /计算其他x值(对于公式中的求和部分,需要调用getmx()函数) result = (float)(bi-getmx(a,x,i,n)/aii); return result; * 解数组Y,计算其中一元素值 * float a 矩阵L * float y 数组Y * int i 数组Y元素序号 * int n 矩阵大小 * return 数组Y的第i个元素值 float gety(float aNN, float bN, float yN, int i, int n) float result; if(i=0) /计算第一个y的值 result = float(bi/aii); else /计算其他y值(对于公式中的求和部分,需要调用getmy()函数) result = float(bi-getmy(a,y,i,n)/aii); return result; int main() float lNN=0; /定义L矩阵 float uNN=0; /定义U矩阵 float yN=0; /定义数组Y float xN=0; /定义数组X float aNN; /定义系数矩阵 float bN; /定义右端项 float sum=0; int i,j,k; int n; int flag=1; while(flag) printf("请输入系数矩阵的大小:"); scanf("%d", &n); if(n>N) printf("矩阵过大!n"); continue; flag=0; printf("请输入系数矩阵值:n"); for(i=0; i<n; i+) for(j=0; j<n; j+) printf("a%d%d: ", i, j); scanf("%f", &aij); printf("请输入右端项数组:n"); for(i=0; i<n; i+) printf("b%d: ", i); scanf("%f", &bi); /*显示原始矩阵*/ printf("n原始矩阵:n"); for(i=0; i<n; i+) for(j=0; j<n; j+) printf("%0.3f ",aij); printf("n"); printf("nn"); /*初始化矩阵l*/ for(i=0; i<n; i+) for(j=0; j<n; j+) if(i=j) lij = 1; /*开始LU分解*/ /*第一步:对矩阵U的首行进行计算*/ for(i=0; i<n; i+) u0i = (float)(a0i/l00); /*第二步:逐步进行LU分解*/ for(i=0; i<n-1; i+) /*对“L列”进行计算*/ for(j=i+1; j<n; j+) for(k=0,sum=0; k<n; k+) if(k != i) sum += ljk*uki; lji = (float)(aji-sum)/uii); /*对“U行”进行计算*/ for(j=i+1; j<n; j+) for(k=0,sum=0; k<n; k+) if(k != i+1) sum += li+1k*ukj; ui+1j = (float)(ai+1j-sum); /*输出矩阵l*/ printf("矩阵L:n"); for(i=0; i<n; i+) for(j=0; j<n; j+) printf("%0.3f ", lij); printf("n"); /*输出矩阵u*/ printf("n矩阵U:n"); for(i=0; i<n; i+) for(j=0; j<n; j+) printf("%0.3f ", uij); printf("n"); /*回代方式计算数组Y*/ for(i=0; i<n; i+) yi = gety(l,b,y,i,n); /*显示数组Y*/ printf("nn数组Y:n"); for(i=0; i<n; i+) printf("y%d = %0.3fn", i+1,yi); /*回代方式计算数组X*/ for(i=n-1; i>=0; i-) xi = getx(u,y,x,i,n); /*显示数组X*/ printf("nn数组X:n"); for(i=0; i<n; i+) printf("x%d = %0.3fn", i+1,xi); return 0;