东北大学数值分析课题2实验报告

数值分析实验报告东北大学一、实验要求1)建立矩阵求逆的算法公式2)编制下三角部分消元，上三角部分消元和对角元单位化等子程序3)应用结构程序设计编程出求n阶矩阵的逆矩阵的通用程序二、实验目的与意义1)通过该课题的实验，掌握求逆矩阵的程序设计方法2)学会利用矩阵求逆运算求解线性方程组的方法，提高数值方法的运用和编程的能力三、实验环境编译环境：Visual C++6.0。

编译语言：C++四、实验过程与分析:(1)求已知矩阵的逆：（以四阶为例）#include<stdio.h>#include <string.h>#define N 4int main(){ float a[N][N];float L[N][N],U[N][N],out[N][N], out1[N][N];float r[N][N],u[N][N];memset( a , 0 , sizeof(a));memset( L , 0 , sizeof(L));memset( U , 0 , sizeof(U));memset( r , 0 , sizeof(r));memset( u , 0 , sizeof(u));int n=N;int k,i,j;int flag=1;float s,t;////////////////////input a matrix////printf("\ninput A=");for(i=0;i<n;i++)for(j=0;j<n;j++)scanf("%f",&a[i][j]);//////////////////figure the input matrix//////////////////////////printf("输入矩阵:\n");for(i=0;i<n;i++){for (j = 0; j < n; j++){printf("%lf ", a[i][j]);}printf("\n");}for(j=0;j<n;j++)a[0][j]=a[0][j]; //计算U矩阵的第一行for(i=1;i<n;i++)a[i][0]=a[i][0]/a[0][0]; //计算L矩阵的第1列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; //计算U矩阵的其他元素}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]; //计算L矩阵的其他元素} }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;}//如果i>j，说明行大于列，计算矩阵的下三角部分，得出L的值，U的//为0else {U[i][j]=a[i][j];if(i==j) L[i][j]=1; //否则如果i<j，说明行小于列，计算矩阵的上三角部分，得出U的//值，L的为0else L[i][j]=0; }} if(U[1][1]*U[2][2]*U[3][3]*U[4][4]==0){flag=0;printf("\n逆矩阵不存在");}if(flag==1){/////////////////////求L和U矩阵的逆for (i=0;i<n;i++) /*求矩阵U的逆*/{u[i][i]=1/U[i][i];//对角元素的值，直接取倒数for (k=i-1;k>=0;k--){s=0;for (j=k+1;j<=i;j++)s=s+U[k][j]*u[j][i];u[k][i]=-s/U[k][k];//迭代计算，按列倒序依次得到每一个值，}}for (i=0;i<n;i++) //求矩阵L的逆{r[i][i]=1; //对角元素的值，直接取倒数，这里为1for (k=i+1;k<n;k++){for (j=i;j<=k-1;j++)r[k][i]=r[k][i]-L[k][j]*r[j][i]; //迭代计算，按列顺序依次得到每一个值}} /////////////////绘制矩阵LU分解后的L和U矩阵///////////////////////printf("\nLU分解后L矩阵:");for(i=0;i<n;i++){ printf("\n");for(j=0;j<n;j++)printf(" %lf",L[i][j]); }printf("\nLU分解后U矩阵:");for(i=0;i<n;i++){ printf("\n");for(j=0;j<n;j++)printf(" %lf",U[i][j]); }printf("\n");////////绘制L和U矩阵的逆矩阵printf("\nL矩阵的逆矩阵:");for(i=0;i<n;i++){ printf("\n");for(j=0;j<n;j++)printf(" %lf",r[i][j]); }printf("\nU矩阵的逆矩阵:");for(i=0;i<n;i++){ printf("\n");for(j=0;j<n;j++)printf(" %lf",u[i][j]); }printf("\n");//验证将L和U相乘，得到原矩阵printf("\nL矩阵和U矩阵乘积\n");for(i=0;i<n;i++){for(j=0;j<n;j++){out[i][j]=0;} }for(i=0;i<n;i++) {for(j=0;j<n;j++){for(k=0;k<n;k++){out[i][j]+=L[i][k]*U[k][j];} } } for(i=0;i<n;i++){for(j=0;j<n;j++){printf("%lf\t",out[i][j]); }printf("\r\n") }//////////将r和u相乘，得到逆矩阵printf("\n原矩阵的逆矩阵:\n");for(i=0;i<n;i++){for(j=0;j<n;j++){out1[i][j]=0;}}for(i=0;i<n;i++){for(j=0;j<n;j++) {for(k=0;k<n;k++){out1[i][j]+=u[i][k]*r[k][j];}} }for(i=0;i<n;i++){for(j=0;j<n;j++){printf("%lf\t",out1[i][j]);}printf("\r\n");}}return 0; }（2）求N阶矩阵的逆：#include<iostream>#include<math.h>using namespace std;class JuZhen{private:double data[10][20];int size;public:void qiuNi();void setSize();void show();void chuShi();};void JuZhen::setSize(){int n;cout<<"请?输?入?矩?阵ó的?阶×:";cin>>n;size=n;}void JuZhen::show(){int i,j;for(i=0;i<size;i++){for(j=size;j<2*size;j++){cout.width(10);cout.flags (ios::right);cout<<data[i][j]<<" ";}cout<<endl;}cout<<"***********************************************"<<endl;cout<<endl;}void JuZhen::chuShi(){int i,j;for(i=0;i<size;i++){cout<<"请?输?入?"<<i+1<<"行D的?元a素?:";for(j=0;j<2*size;j++){if(j<size)cin>>data[i][j];else if(j==i+size)data[i][j]=1.0;elsedata[i][j]=0.0;}}}void JuZhen::qiuNi()int i,j,k;int maxI=0;for(i=1;i<size;i++){if(fabs(data[maxI][0])<fabs(data[i][0])) maxI=i;}if(maxI!=0){double temp;for(j=0;j<2*size;j++){temp=data[0][j];data[0][j]=data[maxI][j];data[maxI][j]=temp;}}double temp2;for(i=0;i<size;i++){if(data[i][i]!=0)temp2=1.0/data[i][i];else{cout<<"此?矩?阵ó无T逆?!"<<endl;return ;}for(j=0;j<2*size;j++)data[i][j]*=temp2;for(j=0;j<size;j++){if(j!=i){double temp3=data[j][i];for(k=0;k<2*size;k++)data[j][k]-=temp3*data[i][k];}}}cout<<"逆?矩?阵ó为a：阰"<<endl;show();}void main()while(1){JuZhen a;a.setSize();a.chuShi();a.qiuNi();}}五、实验结果与分析（1）(2)(3)(4)如图，输入矩阵的阶，然后依次输入矩阵中元素，即可得到逆矩阵，若该矩阵无逆，则输出“此矩阵无逆”。