【連立一次方程式】ガウスの消去法 C++コード付き

今回は連立一次方程式の直接解法の代表である、ガウスの消去法のC++コードを公開します。しくみはおいおい説明するとして、とりあえずコードをおいておきます。反復がi, j, kと出てきて少々ややこしい形になります。

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

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

であり、これを解くと となります。ではコードを以下に示します。

#include <iostream>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <cstdlib>

using namespace std;

int main()
{
	int i,j,k;
	const int n=3;
	double *x=new double[n];
	double *A=new double[n*n];
	double *b=new double[n];
	double m;

	A[0*n+0]=1.;A[0*n+1]=2.;A[0*n+2]=-2.;
	A[1*n+0]=1.;A[1*n+1]=-1.;A[1*n+2]=3.;
	A[2*n+0]=2.;A[2*n+1]=3.;A[2*n+2]=-5.;
	b[0]=3.;b[1]=4.;b[2]=1.;
	
	
	//forward elimination//
	for(i=0;i<n-1;i++)
	{
		for(j=i+1;j<n;j++)
		{
			m=A[j*n+i]/A[i*n+i];
			
			for(k=i+1;k<n;k++)
			{
				A[j*n+k]=A[j*n+k]-m*A[i*n+k];
			}
			
			b[j]=b[j]-m*b[i];
			
			A[j*n+i]=0.;
		}
	}
	
	
	//backward substitution//
	x[n-1]=b[n-1]/A[(n-1)*n+n-1];
	
	for(i=n-2;i>=0;i--)
	{
		for(j=i+1;j<n;j++)
		{
			b[i]-=A[i*n+j]*x[j];
		}
		
		x[i]=b[i]/A[i*n+i];
	}
	
	
	//output//
	ofstream fk;
	fk.open("answer.txt");
	for(i=0;i<n;i++)
	{
		fk<<x[i]<<endl;
		cout<<x[i]<<" ";
	}
	

	delete [] x,A,b;
	return 0;
}