這段代碼前面65行輸出都是對的，就是后面逆矩陣輸出全是0，后面也就不對了
用的lu分解法求逆矩陣#include <stdio.h>
#include<math.h>
int main()
{ int i,j,k,t = 0,s=0,r,n=6; double cond,fanshu1=0,fanshu2=0,l=0,value_l=1,value_u=1; double a[n][n],a_1[n][n]; double L[n][n],U[n][n],L_1[n][n],U_1[n][n]; //生成H矩陣矩陣// for(i=0;i<n;i++) { for(j=0;j<6;j++) { a[i][j]=1.0/(i+j+1); printf("%.10f\n",a[i][j]);//檢查矩陣元素是否正確// } } //矩陣生成// //計算H矩陣的無窮范數// for(i=0;i<n;i++) { for(j=0;j<n;j++) { l=l+a[i][j]; } if(fanshu1<l) fanshu1=l; l=0; } printf("6階H矩陣的范數為：%f\n",fanshu1);//輸出看看h矩陣無窮范數對不對// //計算H矩陣無窮范數結束// //接下來用lu分解求逆矩陣// printf("接下來是l矩陣和u矩陣：\n"); for (i=0;i<n;i++) //計算l矩陣對角線 { L[i][i] = 1.0; } for(k=0;k<n;k++) //計算U// { for (j=k;j<n;j++) { s=0; for (r=0;r<k;r++) { s+=L[k][r]*U[r][k]; } U[k][j]=a[k][j]-s; printf("U[%d][%d]=%f\n",k,j,U[k][j]);//輸出U矩陣看看對不// } for (i=k+1;i<n;i++) { s=0; for(r=0;r<k;r++) { s+=L[i][r]*U[r][k];//按列計算l值 } L[i][k]=(a[i][k]-s)/U[k][k]; printf("L[%d][%d]=%f\n",i,k,L[i][k]);//輸出L矩陣看看對不// } for(i=0;i<n;i++) value_u*=U[i][i]; } //LU分解結束// //下面計算逆矩陣// //l的逆矩陣// for(i=0;i<n;i++) { L_1[i][i]=1; } for(i=1;i<n;i++) { for(j=0;j<i;j++) { t=0; for(k=0;k<i;k++) { t+=L[i][k]*L_1[k][j]; } L_1[i][j]=-t; printf("L_1[%d][%d]=%f\n",i,j,L_1[i][j]); } } //u的逆矩陣// for(i=0;i<n;i++) { U_1[i][i]=1/U[i][i]; } for(i=1;i<n;i++) { for(j=i-1;j>=0;j--) { t=0; for(k=j+1;k<=i;k++) { t+=U[j][k]*U_1[k][i]; } U_1[j][i]=-s/U[j][j]; printf("U_1[%d][%d]=%f\n",i,j,U_1[i][j]); } } //求a的逆矩陣// for(i=0;i<n;i++) { for(j=0;j<n;j++) { t=0; for(r=0;r<n;r++) { t+=U_1[i][r]*L_1[r][j]; } a_1[i][j]=t; } } //求a的逆矩陣的范數// for(i=0;i<n;i++) { for(j=0;j<n;j++) { l=l+a_1[i][j]; printf("%f\n",a_1[i][j]); } if(fanshu2<l) { fanshu2=l; l=0; } printf("fanshu2=%f\n",fanshu2); } cond=fanshu1*fanshu2; printf("cond is：%f\n",cond); }
轉載請註明出處，本文鏈接：https://www.uj5u.com/houduan/85996.html
標籤：C語言