Hatena::ブログ(Diary)

メモ帳Xさん

2011-07-02

LU分解を使ってn元連立1次方程式を解く

LU分解を使ってn元連立1次方程式を解くプログラムをC++で作成した。


LU分解とは、行列を下三角行列L(Lower triangle matrix)と上三角行列U(Upper triangle matrix)の積LUに分解することである。

例えば

A=¥begin{pmatrix}1&2¥¥3&4¥end{pmatrix}

をLU分解すると

L=¥begin{pmatrix}1&0¥¥3&1¥end{pmatrix}
U=¥begin{pmatrix}1&2¥¥0&-2¥end{pmatrix}

となる。(これはあくまで一例であり、この場合A=LUとなるL,Uは無限に存在する)

このLU分解を使ってn元連立1次方程式を解く。


n元連立1次方程式を行列を使って表せば

A¥bf{x}=¥bf{b} (Aはn次の正方行列 ¥bf{x},¥bf{b}はn次の列ベクトル) … (1)

となる。


このときA

A=LU (L,Uはn次の正方行列)

とLU分解する。


すると、(1)の方程式は

LU¥bf{x}=¥bf{b} … (2)

と書き換えられる。


ここで

¥bf{y}=U¥bf{x}

とおけば、(2)より、

L¥bf{y}=¥bf{b}

となる。


このとき、¥bf{L}は下三角行列なので前進代入により、簡単に¥bf{y}が求まる。
さらに、¥bf{U}は上三角行列なので後退代入により、こちらも簡単に¥bf{x}が求まる。

このときの計算量のオーダはO(n^2)である。n元連立1次方程式の解法として有名なガウスの消去法のオーダはO(n^3)なので、LU分解を用いた方が高速である(オーダだけで計算量を完全に比較できるわけではないが、ここでは簡易のためにそうしておく)。

しかし、そのオーダはLU分解に必要なコストを除いたものである。一般に、LU分解をするのに必要な計算量のオーダはガウスの消去法と同じくO(n^3)なので、結果として同じオーダになってしまう。これではLU分解を用いる意味が無いように思える。

もちろん、そのような事は無い。

(1)のような方程式が複数与えられ、そのAがすべて同じ値だったとする。
例えばこんな問題である。

¥begin{pmatrix}1&2¥¥3&4¥end{pmatrix}¥begin{pmatrix}x_1¥¥x_2¥end{pmatrix}=¥begin{pmatrix}1¥¥4¥end{pmatrix}
¥begin{pmatrix}1&2¥¥3&4¥end{pmatrix}¥begin{pmatrix}x_1¥¥x_2¥end{pmatrix}=¥begin{pmatrix}5¥¥6¥end{pmatrix}
¥begin{pmatrix}1&2¥¥3&4¥end{pmatrix}¥begin{pmatrix}x_1¥¥x_2¥end{pmatrix}=¥begin{pmatrix}9¥¥2¥end{pmatrix}

この場合、ガウスの消去法で解くならば全ての問題に対しO(n^3)のオーダの計算量がかかる。
それに対し、LU分解後の行列¥bf{LU}¥bf{A}に対して決まるため、一度計算してしまえば以降同じ¥bf{LU}を使いまわせる。そのため、最初の一問以外はO(n^2)で済む。

これがLU分解を用いる利点である。

LU分解の具体的な方法については省略する(下記のプログラムに書いてはあるが)。

プログラムと問題


プログラム

main.cpp

#include <iostream>

using namespace std;

void LUdecomposition(int dimension, double** A, double** L, double** U){//LU分解
	for(int i=0;i<dimension;i++){
		for(int j=0;j<dimension;j++){
			L[i][j]=0.0;
			U[i][j]=0.0;
			if(i==j)L[i][j]=1.0;
		}
	}
	double sum;
	for(int i=0;i<dimension;i++){
		for(int j=0;j<dimension;j++){
			if(i>j){
				/*L成分を求める*/
				sum=0.0;
				for(int k=0;k<j;k++){
					sum+=L[i][k]*U[k][j];
				}
				L[i][j]=(A[i][j]-sum)/U[j][j];
			}else{
				/*U成分を求める*/
				sum=0.0;
				for(int k=0;k<i;k++){
					sum+=L[i][k]*U[k][j];
				}
				U[i][j]=A[i][j]-sum;
			}
		}
	}
}

void Lforwardsubs(int dimension, double** L, double* b, double* y){//前進代入
	for(int i=0;i<dimension;i++){
		y[i]=b[i];
	}
	for(int i=0;i<dimension;i++){
		y[i]/=L[i][i];
		for(int j=i+1;j<dimension;j++){
			y[j]-=y[i]*L[j][i];
		}
	}
}

void Ubackwardsubs(int dimension, double** U, double* y, double* x){//後退代入
	for(int i=0;i<dimension;i++){
		x[i]=y[i];
	}
	for(int i=dimension-1;i>=0;i--){
		x[i]/=U[i][i];
		for(int j=i-1;j>=0;j--){
			x[j]-=x[i]*U[j][i];
		}
	}
}

int main(){
	int dimension,problem_num;
	double **A,**L,**U,**b,**x,*y;

/*入力を受け取る*/
	cin >> dimension;
	A=new double*[dimension];
	L=new double*[dimension];
	U=new double*[dimension];
	y=new double[dimension];
	for(int i=0;i<dimension;i++){
		A[i]=new double[dimension];
		L[i]=new double[dimension];
		U[i]=new double[dimension];
	}

	for(int i=0;i<dimension;i++){
		for(int j=0;j<dimension;j++){
			cin >> A[i][j];
		}
	}

	cin >> problem_num;
	b=new double*[problem_num];
	x=new double*[problem_num];
	for(int i=0;i<problem_num;i++){
		b[i]=new double[dimension];
		x[i]=new double[dimension];
	}

	for(int i=0;i<problem_num;i++){
		for(int j=0;j<dimension;j++){
			cin >> b[i][j];
		}
	}

/*問題を解く*/
	LUdecomposition(dimension,A,L,U);
	for(int i=0;i<problem_num;i++){
		Lforwardsubs(dimension,L,b[i],y);
		Ubackwardsubs(dimension,U,y,x[i]);
	}

/*解答を出力する*/
	for(int i=0;i<problem_num;i++){
		for(int j=0;j<dimension;j++){
			if(j!=0)cout << " "; 
			cout << x[i][j];
		}
		cout << endl;
	}

/*メモリの解放*/
	for(int i=0;i<dimension;i++){
		delete [] A[i];
		delete [] L[i];
		delete [] U[i];
	}
	for(int i=0;i<problem_num;i++){
		delete [] b[i];
		delete [] x[i];
	}
	delete [] A;
	delete [] L;
	delete [] U;
	delete [] b;
	delete [] x;
	delete [] y;
}

プログラムへの入力は

方程式の次数
係数A(次数の2乗だけ入力)
bの個数
b(上で入力した数だけ入力)

という書式で行う。

下の入力ファイル例を見ればわかるだろう。

入力ファイル(例)

problem.txt

2
2 3
5 7
9
2 3
5 7
11 13
17 19
23 29
31 37
39 41
43 47
53 59


出力ファイル

result.txt

-5 4
-14 11
-38 29
-62 47
-74 57
-106 81
-150 113
-160 121
-194 147


※7/3 少しだけ手直し

スパム対策のためのダミーです。もし見えても何も入力しないでください
ゲスト


画像認証

トラックバック - http://d.hatena.ne.jp/OhXeno/20110702/1309547336