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

数学とか語学とか楽しいよね。ドイツ語とかフランス語とか数値計算とか勉強したことをまとめます。

【数値計算】一次元Poisson方程式(ポアソン方程式)を有限要素法で解く C++コード付き

今回は有限要素法の勉強で最初に解くであろう、一次元Poisson方程式のC++コードを公開します。前に書いた記事「有限要素法のプログラミングを勉強するときにどの順序で学ぶのがよいか?」で紹介した、有限要素法を学ぶ場合に、最初に解くべき方程式です。

普通のGalerkin法で離散化しています。コードは2パターン用意してあります。二次元配列を使うコードと、一次元配列だけで処理しているコードです。二次元配列を使ったほうが行列を足し合わせるときにわかりやすいので、初心者の方はそれでよいと思います。ただ要素数が増えてくると、二次元配列は重いので一次元配列のみで処理するのがおすすめです。

Poisson方程式とは

 \displaystyle \Delta u+f=0

のような二階の偏微分方程式のことを言います。ここで  f は既知の関数です。一次元の場合は

 \displaystyle \frac{ \mathrm{d^2}u }{ \mathrm{d}x^2 }+f=0

となります。以下では領域を  x\in[0,L] とし、  f=1 とします。境界条件は3パターン扱うことができるようにつくってあります。すなわち、

Case 1,  x=0 x=L でDirichlet条件を課す
 u(0) u(L) の値を指定)
Case 2,  x=0 でNeumann条件、 x=L でDirichlet条件を課す
 \frac{ \mathrm{d}u }{ \mathrm{d}x }(0) u(L) の値を指定)
Case 3,  x=0 でDirichlet条件、 x=L でNeumann条件を課す
 u(0) \frac{ \mathrm{d}u }{ \mathrm{d}x }(L) の値を指定)

です。それぞれコード上ではbc=1、bc=2、bc=3に対応しています。連立一次方程式の解法はGauss-Seidel法を用いています。詳しくは以下の記事を参考にして下さい。


まずは二次元配列を使う方からです。

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

using namespace std;

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

//Gauss-Seidel
inline void GS(vector<vector<double> > &A,vector<double> &x,vector<double> &b, int &Node, double &eps)
{
	int i,k;
	double in,sum,max;
	for(k=0;k<100000;k++)
	{
		max=0.0;
		for(i=0;i<Node;i++)
		{
			in=x[i];
			x[i]=(b[i]-(ip(A[i],x,Node)-A[i][i]*x[i]))/A[i][i];
			if(max<abs(in-x[i])) max=abs(in-x[i]);
		}
		if(eps>max) break;
	}
}

inline void mat(vector<vector<double> > &A,vector<double> &b, vector<double> &f, int &Ele, double &dx)
{
	int i;
	for(i=0;i<Ele;i++)
	{
		A[i][i]+=1.0/dx;A[i][i+1]+=-1.0/dx;
		A[i+1][i]+=-1.0/dx;A[i+1][i+1]+=1.0/dx;
		b[i]+=(2.0*f[i]+1.0*f[i+1])*dx/6.0;
		b[i+1]+=(1.0*f[i]+2.0*f[i+1])*dx/6.0;
	}
}

inline void boundary(int &bc,double &D0,double &D1,double &N0,double &N1,vector<vector<double> > &A,vector<double> &b,int &Node)
{
	int i;
	if(bc==1)
	{
		for(i=0;i<Node;i++)
		{
			A[0][i]=0.0;
			A[Node-1][i]=0.0;
		} 
		A[0][0]=1.0;b[0]=D0;
		A[Node-1][Node-1]=1.0;b[Node-1]=D1;
	}
	if(bc==2)
	{
		for(i=0;i<Node;i++)
		{
			A[Node-1][i]=0.0;
		} 
		A[Node-1][Node-1]=1.0;b[Node-1]=D1;
		b[0]+=-N0;
	}
	if(bc==3)
	{
		for(i=0;i<Node;i++)
		{
			A[0][i]=0.0;
		} 
		A[0][0]=1.0;b[0]=D0;
		b[Node-1]+=N1;
	}
}

int main()
{
	int i;
	int Ele=100;
	int Node=Ele+1;
	double L=1.0;
	double dx=L/Ele;
	double eps=pow(2.0,-50);
	int bc=1;//bc=1 both Diriclet bc=2 left Neumann bc=3 right Neumann
	double D0=0.0;
	double D1=2.0; 
	double N0=0.0;
	double N1=2.0;
	vector<double> b(Node,0.);
	vector<double> x(Node,0.);
	vector<double> f(Node,1.);
	vector<vector<double> > A(Node,(vector<double>(Node,0.)));
	
	ofstream fk;
	fk.open("answer.txt");
	
	mat(A,b,f,Ele,dx);
	boundary(bc,D0,D1,N0,N1,A,b,Node);
	GS(A,x,b,Node,eps);
	
	for(i=0;i<Node;i++) cout<<x[i]<<endl;
	for(i=0;i<Node;i++) fk<<dx*i<<" "<<x[i]<<endl;

	return 0;
}


次に一次元配列だけで処理している方です。今回は行列が三重対角行列になるので、vectorのDに対角要素を、Lに対角要素の左隣の要素を、Rに対角要素の右隣の要素を格納しています。

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

using namespace std;

//Gauss-Seidel
inline void GS(vector<double> &D,vector<double> &L,vector<double> &R,vector<double> &x,vector<double> &b, int &Node, double &eps)
{
	int i,k;
	double in,sum,max;
	for(k=0;k<100000;k++)
	{
		max=0.0;
		in=x[0];
		x[0]=(b[0]-R[0]*x[1])/D[0];
		if(max<abs(in-x[0])) max=abs(in-x[0]);
		for(i=1;i<Node-1;i++)
		{
			in=x[i];
			x[i]=(b[i]-(L[i]*x[i-1]+R[i]*x[i+1]))/D[i];
			if(max<abs(in-x[i])) max=abs(in-x[i]);
		}
		in=x[Node-1];
		x[Node-1]=(b[Node-1]-L[Node-1]*x[Node-2])/D[Node-1];
		if(max<abs(in-x[Node-1])) max=abs(in-x[Node-1]);
		if(eps>max) break;
	}
}

inline void mat(vector<double> &D,vector<double> &L,vector<double> &R,vector<double> &b, vector<double> &f, int &Ele, double &dx)
{
	int i;
	for(i=0;i<Ele;i++)
	{
		D[i]+=1.0/dx;R[i]+=-1.0/dx;
		L[i+1]+=-1.0/dx;D[i+1]+=1.0/dx;
		b[i]+=(2.0*f[i]+1.0*f[i+1])*dx/6.0;
		b[i+1]+=(1.0*f[i]+2.0*f[i+1])*dx/6.0;
	}
}

inline void boundary(int &bc,double &D0,double &D1,double &N0,double &N1,vector<double> &D,vector<double> &L,vector<double> &R,vector<double> &b,int &Node)
{
	if(bc==1)
	{
		D[0]=1.0;R[0]=0.0;b[0]=D0;
		L[Node-1]=0.0;D[Node-1]=1.0;b[Node-1]=D1;
	}
	if(bc==2)
	{
		L[Node-1]=0.0;D[Node-1]=1.0;b[Node-1]=D1;
		b[0]+=-N0;
	}
	if(bc==3)
	{
		D[0]=1.0;R[0]=0.0;b[0]=D0;
		b[Node-1]+=N1;
	}
}

int main()
{
	int i;
	int Ele=100;
	int Node=Ele+1;
	double LL=1.0;
	double dx=LL/Ele;
	double eps=pow(2.0,-50);
	int bc=1;//bc=1 both Diriclet bc=2 left Neumann bc=3 right Neumann
	double D0=0.0;
	double D1=2.0; 
	double N0=0.0;
	double N1=2.0;
	vector<double> b(Node,0.);
	vector<double> x(Node,0.);
	vector<double> f(Node,1.);
	vector<double> D(Node,0.);
	vector<double> L(Node,0.);
	vector<double> R(Node,0.);
	
	ofstream fk;
	fk.open("answer.txt");
	
	mat(D,L,R,b,f,Ele,dx);
	boundary(bc,D0,D1,N0,N1,D,L,R,b,Node);
	GS(D,L,R,x,b,Node,eps);
	
	for(i=0;i<Node;i++) cout<<x[i]<<endl;
	for(i=0;i<Node;i++) fk<<dx*i<<" "<<x[i]<<endl;

	return 0;
}

コードの流れですが、まず要素の分割数と境界条件の型(bcの値による)を決め、境界条件の設定(値をどうするか)を行い、関数 mat により行列を作成、関数 boundary で境界条件を行列に反映させ、関数 GS(Gauss-Seidel法)で作成した行列を解き、出力しています。出力はテキストファイル "answer.txt" に出てきます。左の行が  x 座標で右の行が求めている関数値  u です。

今回はコードがメインで離散化の詳細を示していませんが、必ず書くのでお待ちください。