LU分解法 逆行列の計算
Posted: 2009年1月12日(月) 02:00
また同じような質問なので、すみません。
LU分解を用いて、3次正方行列の逆行列の計算をする問題なのですが、分かりませんでした。
まず、上手くLとUに分解できたのですか、その後のLの逆行列とUの逆行列を出すことができませんでした。
そして、問題の答えである、Lの逆行列とUの逆行列の掛算がうまくできませんでした。
LとUの逆行列ができないので、それはそうなのですが、すみません。
LとUの分解を以下に載せておきます。
同じような問題である、3次正方行列の方程式を解く問題は解けましたので、以下に載せておきます。
非常に長くなりましたが、暇のある人は挑戦してみてください。お願いします。
LU分解を用いて、3次正方行列の逆行列の計算をする問題なのですが、分かりませんでした。
まず、上手くLとUに分解できたのですか、その後のLの逆行列とUの逆行列を出すことができませんでした。
そして、問題の答えである、Lの逆行列とUの逆行列の掛算がうまくできませんでした。
LとUの逆行列ができないので、それはそうなのですが、すみません。
LとUの分解を以下に載せておきます。
#include <stdio.h>
#include <math.h>
#define N 3
float A[N][N]={{2.0,1.0,1.0},{1.0,5.0,-1.0},{-1.0,1.0,-3.0}};
int main(void){
int i,j,k,p,q;
float L[N][N]={0.0},U[N][N]={0.0};
float lu;
for(j=0;j<N;j++) L[j][j]=1.0;
for(j=0;j<N;j++) U[0][j]=A[0][j];
for(j=1;j<N;j++) L[j][0]=A[j][0]/U[0][0];
for(j=1;j<N-1;j++){
for(q=j;q<N;q++){
lu=0.0;
for(k=0;k<j;k++) lu+=L[j][k]*U[k][q];
U[j][q]=A[j][q]-lu;
}
for(p=j+1;p<N;p++){
lu=0.0;
for(k=0;k<j;k++) lu+=L[p][k]*U[k][j];
L[p][j]=(A[p][j]-lu)/U[j][j];
}
}
lu=0.0;
for(k=0;k<N-1;k++)lu+=L[N-1][k]*U[k][N-1];
U[N-1][N-1]=A[N-1][N-1]-lu;
for(i=0;i<N;i++){
printf("[ ");
for(j=0;j<N;j++)printf("%9f ",A[j]);
printf("]\n");
}
printf("のLU分解は\n");
for(i=0;i<N;i++){
printf("[ ");
for(j=0;j<N;j++)printf("%9f ",L[j]);
printf("]\n");
}
printf("\n");
for(i=0;i<N;i++){
printf("[ ");
for(j=0;j<N;j++)printf("%9f ",U[j]);
printf("]\n");
}
return(0);
}同じような問題である、3次正方行列の方程式を解く問題は解けましたので、以下に載せておきます。
#include <stdio.h>
#include <math.h>
#define N 3
float A[N][N]={{2.0,1.0,1.0},{1.0,5.0,-1.0},{-1.0,1.0,-3.0}};
float B[N]={3.0,-6.0,-8.0};
int main(void){
int j,k,p,q;
float L[N][N]={0.0},U[N][N]={0.0};
float lu;
float x[N],y[N];
for(j=0;j<N;j++) L[j][j]=1.0;
for(j=0;j<N;j++) U[0][j]=A[0][j];
for(j=1;j<N;j++) L[j][0]=A[j][0]/U[0][0];
for(j=1;j<N-1;j++){
for(q=j;q<N;q++){
lu=0.0;
for(k=0;k<j;k++) lu+=L[j][k]*U[k][q];
U[j][q]=A[j][q]-lu;
}
for(p=j+1;p<N;p++){
lu=0.0;
for(k=0;k<j;k++) lu+=L[p][k]*U[k][j];
L[p][j]=(A[p][j]-lu)/U[j][j];
}
}
lu=0.0;
for(k=0;k<N-1;k++)lu+=L[N-1][k]*U[k][N-1];
U[N-1][N-1]=A[N-1][N-1]-lu;
y[0]=B[0];
for(k=1;k<N;k++){
lu=0.0;
for(j=0;j<k;j++)lu+=L[k][j]*y[j];
y[k]=B[k]-lu;
}
x[N-1]=y[N-1]/U[N-1][N-1];
for(k=N-2;k>=0;k--){
lu=0.0;
for(j=N-1;j>k;j--)lu+=U[k][j]*x[j];
x[k]=(y[k]-lu)/U[k][k];
}
printf("解:x=%f y=%f z=%f\n",x[0],x[1],x[2]);
return(0);
}非常に長くなりましたが、暇のある人は挑戦してみてください。お願いします。