数学とか語学とか楽しいよね

フランス語、ドイツ語、ロシア語、アラビア語、オランダ語、英語、スペイン語、ラテン語とか数学とか数値計算(有限要素法、有限体積法、差分法、格子ボルツマン法、数理最適化、C++コード付き)とか勉強したことをまとめます。右のカテゴリーから興味のある記事を探してください。最近はクラシックの名演も紹介しています。Amazonアソシエイトを使用しています。

【数値計算】LU分解(Doolittle法)C++コード付き

今回は連立一次方程式に対する直接法である、LU分解(Doolittle法)のC++コードを公開します。

今回解く連立一次方程式は

 Ax=b

とあらわされます。ここで、  A 4 \times 4 の行列、 x b 4 次元のベクトルです。 A b が既知で、 x が未知です。具体的には

 \left( \begin{array}{ccc} 8 & 16 & 24 & 32 \\ 2 & 7 & 12 & 17 \\ 6 & 17 & 32 & 59 \\ 7 & 22 & 46 & 105 \end{array} \right) \left( \begin{array}{c} x_1 \\ x_2 \\ x_3 \\ x_4 \end{array} \right) = \left( \begin{array}{c} 80 \\ 38 \\ 114 \\ 180 \end{array} \right)

であり、これを解くと  x_1=1, x_2=1, x_3=1, x_4=1 となります。プログラムでは  A b はテキストファイルから読み込ませています。また、LU分解の結果は一つの配列(行列)LUに格納しています。ではコードを以下に示します。

#include <iostream>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <stdio.h>
#include <time.h>
#include <stdlib.h> 

using namespace std;

//matix vector multiplication Ax=b
inline void mvm(double A[],double x[],double b[], int size)
{
	int i,j;
	
	for(i=0;i<size;i++)
	{
		b[i]=0.;
		for(j=0;j<size;j++)
		{
			b[i]+=A[i*size+j]*x[j];
		}
	}
}

//inner product//
inline double ip(double a[], double b[], int size)
{
	int i;
	
	double value=0.;
	for(i=0;i<size;i++) value+=a[i]*b[i];
	return value;
}

//input//
inline void input(double A[], double b[], int size)
{
	int i,j;
	
	ifstream fo;
	fo.open("A.txt");
	for(i=0;i<size;i++)
	{
		for(j=0;j<size;j++)
		{
			fo>>A[i*size+j];
		}
	}
	fo.close();
	
	fo.open("b.txt");
	for(i=0;i<size;i++)
	{
		fo>>b[i];
	}
	fo.close();
}

//matrix output//
inline void matout(double A[], int size)
{
	int i,j;
	
	for(i=0;i<size;i++)
	{
		for(j=0;j<size;j++)
		{
			cout<<A[i*size+j]<<" ";
		}
		cout<<endl;
	}
	cout<<endl;
}

//vector output//
inline void vecout(double b[], int size)
{
	int i;
	
	for(i=0;i<size;i++)
	{
		cout<<b[i]<<endl;
	}
	cout<<endl;
}

//LU decomposition//
inline void LUd(double A[], double LU[], int size)
{
	int i,j,k;
	double sum;
	
	for(i=0;i<size;i++)
	{
		for(j=0;j<size;j++)
		{
			if(i>j)//l_ij
			{
				sum=0.0;
				for(k=0;k<j;k++) sum+=LU[i*size+k]*LU[k*size+j];
				LU[i*size+j]=(A[i*size+j]-sum)/LU[j*size+j];
			}
			else//u_ij
			{
				sum=0.0;
				for(k=0;k<i;k++) sum+=LU[i*size+k]*LU[k*size+j];
				LU[i*size+j]=A[i*size+j]-sum;
			}
		}
	}
}

//forward substitution and backward substitution//
void inline fb(double LU[], double x[], double y[], double b[], int size)
{
	int i,j;
	double sum;
	
	y[0]=b[0];
	
	//forward//
	for(i=1;i<size;i++)
	{
		sum=0.0;
		for(j=0;j<i;j++)
		{
			sum+=LU[i*size+j]*y[j];
		}
		y[i]=b[i]-sum;
	}
	
	//backward//
	x[size-1]=y[size-1]/LU[(size-1)*size+size-1];
	
	for(i=size-2;i>=0;i--)
	{
		sum=1.0*x[i];
		for(j=i+1;j<size;j++)
		{
			sum+=LU[i*size+j]*x[j];
		}
		x[i]=(y[i]-sum)/LU[i*size+i];
	}
}

main()
{
	int size=4;
	
	double *b=new double[size];
	double *x=new double[size];
	double *A=new double[size*size];
	double *LU=new double[size*size];
	double *y=new double [size];
	
	//input//
	input(A,b,size);
	
	//LU decomposition//
	LUd(A,LU,size);
	
	//forward substitution and backward substitution//
	fb(LU,x,y,b,size);
	
	//output//
	matout(LU,size);
	vecout(x,size);
	
	delete [] b,x,A,LU,y;

    return 0;
}


参考