LU分解を用いた逆行列の求め方
お世話になっております。
C++でLU分解と連立一次方程式の解を求める関数SLE_bu_LUを利用して、逆行列の求める関数InvMatを作成中です。
「行列Bのj列目の列ベクトルをbj とすれば，行列Cのj列目の列ベクトル
xj は連立１次方程式 Axj = bj の解として求められる．これが各列について成立するので，
求めたx1, x2, x3，････, xmを各列としてC=[x1 x2 x3 ････ xm] を計算することがで
きる。その際，係数行列Aが同一の連立１次方程式をm 組解いて行列C が得られること
がわかる。従って，初めに１回だけAを LU 分解しておいて，前進代入と後退代入のプロセ
スをm回行うことによって行列C が求められる」
と議題には書かれており、Cを求める逆行列、Bを単位行列としてプログラミングを書いていたのですが、
関数InvMatではCとBのm列成分を関数SLE by LUに引数として渡すという方針であっていますでしょうか？
（そうなるとm回関数SLE_by_LUを呼び出し、m回LU分解してしまうので間違っていると思いました...）
関数SLE_bu_LUは正しく動き、追加した関数は一番下のInvMatのみでここが怪しいと思います。
よろしくお願いします。
#include <iostream>
#include <iomanip>
using namespace std;
const int N = 4; /* 連立方程式の元数 */
void SLE_by_LU(double *,double [][N],double *); // LU 分解により連立１次方程式
// を解く関数のプロトタイプ宣言
void InvMat(double [][N],double [][N]); //逆行列を求める関数のプロトタイプ宣言
int main()
{
double x[N];
/* 入力データ （演習問題 2.1 の場合） */
double a[N][N]={{1.0,2.0,-12.0,8.0},{5.0,4.0,7.0,-2.0},
{-3.0,7.0,9.0,5.0},{6.0,-12.0,-9.0,3.0}};
double ai[N][N]={}; //逆行列用
double b[N]={27.0, 4.0,11.0,49.0};
SLE_by_LU(x,a,b); // LU 分解により連立１次方程式を解く関数の呼び出し
/* データ出力 */
for(int i=0; i<N; i++) {cout << 'x' << i+1 << '=' << setw(12) <<
fixed << setprecision(9) << x[i] << endl;
}
InvMat(a,ai);//逆行列の計算関数呼び出し
return 0;
}
/****** LU 分解により連立１次方程式を解く関数プログラム **********/
/* 入力 ： a[N][N]=係数行列，b[N]=右辺ベクトル，N=元数 */
/* 出力 ： x[N]=方程式の解 */
/******************************************************************/
void SLE_by_LU(double x[],double a[][N],double b[])
{ int i, j, k;
double l[N][N], u[N][N], y[N]; // LU 分解用配列
double sigma; // ∑演算用作業変数
/* LU 分解 */
for(j=0; j<N; j++) u[0][j]=a[0][j]; // ステップ 1.1
for(j=1; j<N; j++) l[j][0]=a[j][0]/u[0][0]; // ステップ 1.2
for(k=1; k<N; k++) // ステップ 2 の反復
{
for(j=k; j<N; j++) // ステップ 2.1
{
u[k][j]=a[k][j];
for(i=0; i<k; i++) u[k][j]-=l[k][i]*u[i][j];
}
for(j=k+1; j<N; j++) // ステップ 2.2
{
sigma=a[j][k];
for(i=0; i<k; i++) sigma-=l[j][i]*u[i][k];
l[j][k]=sigma/u[k][k];
}
}
/* 前進代入 */
for(i=0; i<N; i++) // i=1,2,･･･,n
{
y[i]=b[i]; // yi ← bi
for(j=0; j<i; j++) y[i]-=l[i][j]*y[j]; // yi=bi-∑lijyj
}
/* 後退代入 */
for(k=N-1; k>=0; k--)
{
sigma=y[k];
for(j=k+1; j<N; j++) sigma-=u[k][j]*x[j];
x[k]=sigma/u[k][k];
}
}
void InvMat(double a[][N], double ai[][N]) //A 配列a[][N]に格納 逆行列AI 配列ai[][N]に格納
{//AC=BでBが単位行列Eの時、Cは逆行列AIとなる
double array[N][N];//単位行列格納
double array_m[N];//単位行列のm列目格納
double ai_m[N];//AIのm列目格納
for(int i = 0; i <N; i++){//単位行列を作る
for(int j = 0; j <N; j++){
if(i == j )
{array[i][j] = 1.0;}
else
{array[i][j] = 0.0;}
}
}
//以下A.14より
for(int i=0;i<N;i++){//m列目処理
for(int j=0;j<N;j++){
array_m[i]=array[j][i];//m番目の単位行列の列成分
ai_m[i]=ai[j][i];//m番目のAI行列の列成分
SLE_by_LU(ai_m,a,array_m);
cout<<ai[i][j];
}
}
}
0
