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

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

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

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


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

 Ax=b

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

 \left( \begin{array}{ccc} 1 & 2 & -2 \\ 1 & -1 & 3 \\ 2 & 3 & -5 \end{array} \right) \left( \begin{array}{c} x_1 \\ x_2 \\ x_3 \end{array} \right) = \left( \begin{array}{c} 3 \\ 4 \\ 1 \end{array} \right)

であり、これを解くと  x_1=1, x_2=3, x_3=2 となります。ではコードを以下に示します。

#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;
}