实验六报告6LU分解

实验六解线性方程组的直接法
——LU分解
实验目的:
1、了解LU分解的数学依据与基本原理；
2、能熟练地应用、构造矩阵的LU分解程序求解方程组。

实验内容:
根据LU分解的原理与公式自己填写
实验题目：
（1）编写程序写出下面系数矩阵的LU分解
（2）编程序，分别用LU分解法与主元素消去法解方程组
第一题：
#include<stdio.h>
#include<stdlib.h>
void main()
{
float a[3][3]={{2,2,3},{4,5,7},{-2,4,5}};
float L[3][3],U[3][3];
int n=3;
int k,i,j;
float s,t;
for(j=0;j<n;j++)
a[0][j]=a[0][j];
for(i=1;i<n;i++)
a[i][0]=a[i][0]/a[0][0];
for(k=1;k<n;k++)
{
for(j=k;j<n;j++)
{
s=0;
for (i=0;i<k;i++)
s=s+a[k][i]*a[i][j];
a[k][j]=a[k][j]-s;
}
for(i=k+1;i<n;i++)
{
t=0;
for(j=0;j<k;j++)
t=t+a[i][j]*a[j][k];
a[i][k]=(a[i][k]-t)/a[k][k];
}
}
for(i=0;i<n;i++)
for(j=0;j<n;j++)
{
if(i>j)
{
L[i][j]=a[i][j]; U[i][j]=0;
}
else
{
U[i][j]=a[i][j];
if(i==j)
L[i][j]=1;
else
L[i][j]=0;
}
}
printf("L[3][3]=\n");
for(i=0;i<n;i++)
for(j=0;j<n;j++)
printf(" %0.0f",L[i][j]);
printf("\n");
}
printf("U[3][3]=\n");
for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
printf(" %0.0f",U[i][j]);
printf("\n");
}
}
第二题：
#include <stdio.h>
#include <stdlib.h>
#define N 10 //矩阵大小范围
//求x回调函数
float getmx(float a[N][N], float x[N], int i, int n) {
float mx = 0;
int r;
for(r=i+1; r<n; r++)
{
mx += a[i][r] * x[r];
}
return mx;
}
//求y回调函数
float getmy(float a[N][N], float y[N], int i, int n) {
float my = 0;
int r;
for(r=0; r<n; r++)
{
if(i != r) my += a[i][r] * y[r];
}
return my;
}
//计算x
float getx(float a[N][N], float b[N], float x[N], int i, int n) {
float result;
if(i==n-1) //计算最后一个x的值
result = (float)(b[i]/a[n-1][n-1]);
else
result = (float)((b[i]-getmx(a,x,i,n))/a[i][i]);
return result;
}
//计算y
float gety(float a[N][N], float b[N], float y[N], int i, int n) {
float result;
if(i==0) //计算第一个y的值
result = float(b[i]/a[i][i]);
else
result = float((b[i]-getmy(a,y,i,n))/a[i][i]);
return result;
}
void main()
{ float l[N][N]={0}; //定义L矩阵
float u[N][N]={0}; //定义U矩阵
float y[N]={0}; //定义数组Y
float x[N]={0}; //定义数组X
float a[N][N]; //定义系数矩阵
float b[N]; //定义右端项
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=1; i<n+1; i++)
{
for(j=1; j<n+1; j++)
{
printf("a[%d][%d]: ", i, j);
scanf("%f", &a[i][j]);
}
}
printf("请输入右端项数组：\n");
for(i=1; i<n+1; i++)
{
printf("b[%d]: ", i);
scanf("%f", &b[i]);
}
//初始化矩阵l
for(i=0; i<n; i++)
{
for(j=0; j<n; j++)
{
if(i==j) l[i][j] = 1;
}
}
//LU分解
for(i=0; i<n; i++)
{
u[0][i] = (float)(a[0][i]/l[0][0]);
}
for(i=0; i<n-1; i++)
{
for(j=i+1; j<n; j++)
{
for(k=0,sum=0; k<n; k++)
{
if(k != i) sum += l[j][k]*u[k][i]; }
l[j][i] = (float)((a[j][i]-sum)/u[i][i]); }
for(j=i+1; j<n; j++)
{
for(k=0,sum=0; k<n; k++)
{
if(k != i+1) sum += l[i+1][k]*u[k][j]; }
u[i+1][j] = (float)((a[i+1][j]-sum));
}
}
//回代方式计算数组
for(i=0; i<n; i++)
{
y[i] = gety(l,b,y,i,n);
}
//回代方式计算数组X
for(i=n-1; i>=0; i--)
{
x[i] = getx(u,y,x,i,n);
}
//显示数组X即结果
printf("\n\n方程组解X：\n");
for(i=0; i<n; i++)
{
printf("x%d = %0.3f\n", i+1,x[i]);
}
}。