LU分解法

学院名称 学生姓名 课程名称 计算机科学与技术 王胜涛 数值计算方法 专业班级 实验日期 实验题目 计算机0901班 2011-6-13 LU分解法 学号 成绩 2009001448 一、实验目的及任务 线性方程组的直接解法 合理利用Gauss消元法求解下列方程组： 484125x12x772102 836x3761120x4321 二、源代码 #include \"stdafx.h\" #include \"stdio.h\" #include \"stdlib.h\" #define withresult 1 #define withoutresult 0 double A [ 15 ] [ 15 ] , B[ 15 ] , R[15] , Y[ 15 ] ; 所用的存储空间 int n ; //存储系数矩阵大小的变量 void DisplayA( unsigned style ) 打印线性方程组的函数 { printf( \"\\n\" ) ; for ( int i = 0 ; i < n ; i ++ ) 印线性方程组 { for ( int j = 0 ; j < n ; j ++ ) { printf( \"a[%d][%d] = %6.2lf \" , i , j , A [ i ] [ j ] ) ; } printf ( \"b[%d]=%6.2lf\\n\" , i , B [ i ] )

//存储线性方程组 // //打 ; } //打印解向量 if ( style == withresult ) { printf (\"\\n解向量为：\\n\") ; for ( int i = 0 ; i < n ; i ++ ) { printf ( \"R[%d] = %6.2lf \\n\" , i , R [ i ] ) ; } } printf ( \"\\n\" ) ; } int _tmain(int argc, _TCHAR* argv[]) { //从文件中读入线性方程组 FILE * input ; input = fopen (\"LU.txt.txt\" ,\"r\" ) ; //打开存储有线性方程组的文件 //读入矩阵大小 fscanf( input , \"%d\" , & n ) ; //读入系数矩阵 for ( int i = 0 ; i < n ; i ++ ) //打印线性方程组 { for ( int j = 0 ; j < n ; j ++ ) { fscanf( input ,\"%lf\" , & A [ i ] [ j ] ) ; } fscanf ( input , \"%lf\" , & B [ i ] ) ; } DisplayA ( withoutresult) ; //下面开始LU分解 for ( int k=0 ; k < n ; k ++ ) { if ( !k ) continue ; //第一行的值没有改变 else { for( int i = 0 ; i < k ; i++ ) {

if ( i == 0 ) //计算第一列的值 { A[ k ] [ i ] = A[ k ] [ i ] / A [ 0 ] [ 0 ] ; } else { double sum = 0 ; for( int m = 0 ; m < i ; m ++ ) { sum += A [ k ] [ m ] * A [ m ] [ i ] ; } A [ k ] [ i ] = ( A [ k ] [ i ] - sum ) / A [ i ][ i ] ; } } for ( int i = k ; i < n ; i ++ ) { double sum = 0 ; for ( int m = 0 ; m < k ; m ++) { sum += A[ k ] [ m ] * A [ m ] [ i ]; } A[ k ] [ i ] = A [ k ] [ i ] - sum ; } } } //LU分解完毕，开始求解答案 //先求解Y解向量 for ( int i = 0 ; i < n ; i ++ ) { if( ! i ) Y[ i ] = B [ i ] ; else { double sum = 0 ; for( int j = 0 ; j < i ; j ++ ) { sum += A [ i ] [ j ] * Y[ j ] ; } Y[ i ] = B [ i ] - sum ; } } //计算最后的结果解向量 for ( int i = n - 1 ; i >= 0 ; i -- ) {

if ( i == n -1 ) { R [ i ]= Y [ i ] / A[ i ] [ i ] ; } else { double sum = 0 ; for ( int j = i + 1 ; j < n ; j ++ ) { sum += A [ i ] [ j ] * R [ j ] } R [ i ] = ( Y [ i ] - sum ) / A [ i ] [ i ] } } DisplayA ( withresult ) ; system ( \"pause\" ) ; return 0 ; } 三、实试验结果及分析 ; ; 分析： 实验地点

综合楼六层606室 指导教师 王峥