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

ドイツ語、オランダ語、英語、フランス語、ラテン語とか数値計算とか勉強したことをまとめます。右のカテゴリーから興味のある記事を探してください。最近はクラシックの名演も紹介しています。

【数値計算】一次元非定常拡散方程式を有限要素法で解く C++コード付き

今回は一次元非定常拡散方程式を有限要素法で解いていきます。C++コード付きです。一次元非定常拡散方程式は、前に書いた記事「有限要素法のプログラミングを勉強するときにどの順序で学ぶのがよいか?」で紹介した、有限要素法を学ぶ場合に三番目に解くべき方程式です。

一次元非定常拡散方程式とは

 \displaystyle \frac{\partial u}{\partial t}=D\frac{\partial^2 u}{\partial x^2} in  [0, L]

のような偏微分方程式のことをいいます。ここで、 u(t,x) は熱や溶質の濃度、 D は拡散の効果を表す拡散係数と解釈することができます。熱またはある溶質が時間が経過するにつれてどのような分布をとるか(拡散していくか)求めるのがこの問題です。 u は時間方向と空間方向ともに変化するので二変数関数  u(t,x) となっています。前に説明した一次元Poisson方程式を時間発展させたものにあたります。今回も二次元配列を使うコードと一次元配列だけで処理するコードの両方を公開します。境界条件は前回と同じく3パターン用意しました。すなわち、

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

です。それぞれコード上ではbc=1、bc=2、bc=3に対応しています。連立一次方程式の解法はGauss-Seidel法を用いています。空間方向の離散化はGalerkin法、時間方向は  \theta 法で離散化しており、今回は  \theta=0.5 のCrank-Nicolson法を用いています。ちなみに、 \theta=0 のとき前進Euler法(陽解法、explicit scheme)に、 \theta=1 のときに後退Euler法(完全陰解法、implicit scheme)になります。これらについても今後まとめる予定です。

計算条件については以下のように設定しています。まず、拡散方程式を解く領域は  x\in[0,1]、要素数は100、節点数は101で等分割( \Delta x=1/100=0.01)、時間刻みは  \Delta t=0.001、計算終了時刻は  t=10、拡散係数は  D=0.01 \theta=0.5境界条件はbc=1(両側Dirichlet条件)で u(t,0)=0 u(t,1)=0、初期条件は頂点の値が0.5の三角形(下の図参照)です。初期時刻に三角形のような熱の分布を与え、棒の両側をずっと0℃で冷やしたときに、棒の熱はある時間位どのような分布をとるか、という問題です。答えは"answer.txt"に出力されます。横方向が  x に対応しています。一行目が  x 座標、二行目が初期条件(時刻  t=0 における  u の分布、つまり  u(0,x))、三行目が時刻  t=1 における  u の分布( u(1,x))、…となっています。


まずは二次元配列を使う方です。
(準備中)


次に一次元配列だけで処理する方です。

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

using namespace std;

//Gauss-Seidel
inline void GS(vector<double> &AD,vector<double> &AL,vector<double> &AR,vector<double> &BD,vector<double> &BL,vector<double> &BR,vector<double> &x,vector<double> &xx,vector<double> &b, int &Node, double &eps, int &bc, double &D0,double &D1,double &N0,double &N1)
{
	int i,k;
	double in,max;
	for(k=0;k<100000;k++)
	{
		if(bc==1)
		{
			b[0]=D0;
			b[Node-1]=D1;
		}
		if(bc==2)
		{
			b[Node-1]=D1;
			b[0]=BD[0]*x[0]+BR[0]*x[1]-N0;
		}
		if(bc==3)
		{
			b[0]=D0;
			b[Node-1]=BL[Node-1]*x[Node-2]+BD[Node-1]*x[Node-1]+N1;
		}
		
		for(i=1;i<Node-1;i++)
		{
			b[i]=BL[i]*x[i-1]+BD[i]*x[i]+BR[i]*x[i+1];
		}
		
		max=0.0;
		in=xx[0];
		xx[0]=(b[0]-AR[0]*xx[1])/AD[0];
		if(max<abs(in-xx[0])) max=abs(in-xx[0]);
		for(i=1;i<Node-1;i++)
		{
			in=xx[i];
			xx[i]=(b[i]-(AL[i]*xx[i-1]+AR[i]*xx[i+1]))/AD[i];
			if(max<abs(in-xx[i])) max=abs(in-xx[i]);
		}
		in=xx[Node-1];
		xx[Node-1]=(b[Node-1]-AL[Node-1]*xx[Node-2])/AD[Node-1];
		if(max<abs(in-xx[Node-1])) max=abs(in-xx[Node-1]);
		if(eps>max) break;
	}
}

inline void mat(vector<double> &AD,vector<double> &AL,vector<double> &AR,vector<double> &BD,vector<double> &BL,vector<double> &BR,vector<double> &b, vector<double> &f, int &Ele, double &dx, double &Diff, double &dt, double &theta)
{
	int i;
	for(i=0;i<Ele;i++)
	{
		//A//
		AD[i]+=dx*2.0/6.0/dt+theta*Diff/dx;
		AR[i]+=dx*1.0/6.0/dt+theta*-Diff/dx;
		AL[i+1]+=dx*1.0/6.0/dt+theta*-Diff/dx;
		AD[i+1]+=dx*2.0/6.0/dt+theta*Diff/dx;
		
		//B//
		BD[i]+=dx*2.0/6.0/dt-(1.0-theta)*Diff/dx;
		BR[i]+=dx*1.0/6.0/dt-(1.0-theta)*-Diff/dx;
		BL[i+1]+=dx*1.0/6.0/dt-(1.0-theta)*-Diff/dx;
		BD[i+1]+=dx*2.0/6.0/dt-(1.0-theta)*Diff/dx;
	}
}

inline void boundary(int &bc,vector<double> &AD,vector<double> &AL,vector<double> &AR,vector<double> &BD,vector<double> &BL,vector<double> &BR,vector<double> &b,int &Node)
{
	if(bc==1)
	{
		AD[0]=1.0;AR[0]=0.0;BD[0]=1.0;BR[0]=0.0;
		AL[Node-1]=0.0;AD[Node-1]=1.0;BL[Node-1]=0.0;BD[Node-1]=1.0;
	}
	if(bc==2)
	{
		AL[Node-1]=0.0;AD[Node-1]=1.0;BL[Node-1]=0.0;BD[Node-1]=1.0;
	}
	if(bc==3)
	{
		AD[0]=1.0;AR[0]=0.0;BD[0]=1.0;BR[0]=0.0;
	}
}

int main()
{
	int i,j;
	int Ele=100;
	int Node=Ele+1;
	double LL=1.0;
	double dx=LL/Ele;
	double dt=0.001;
	double NT=10000;
	double eps=pow(2.0,-50);
	double Diff=0.01;
	double theta=0.5;
	int bc=1;//bc=1 both Diriclet bc=2 left Neumann bc=3 right Neumann
	double D0=0.0;
	double D1=0.0;
	double N0=0.0;
	double N1=0.0;
	vector<double> b(Node,0.);
	vector<double> x(Node,0.);
	vector<double> xx(Node,0.);
	vector<double> f(Node,0.);
	vector<double> AD(Node,0.);
	vector<double> AL(Node,0.);
	vector<double> AR(Node,0.);
	vector<double> BD(Node,0.);
	vector<double> BL(Node,0.);
	vector<double> BR(Node,0.);
	
	ofstream fk;
	fk.open("answer.txt");
	
	for(i=0;i<Node;i++) fk<<dx*float(i)<<" ";
	fk<<endl;
	
	//initial condition//
	for(i=0;i<Node;i++)
	{
		x[i]=min(dx*float(i),1.0-dx*float(i));
		fk<<x[i]<<" ";
	}
	fk<<endl;
	
	
	mat(AD,AL,AR,BD,BL,BR,b,f,Ele,dx,Diff,dt,theta);
	boundary(bc,AD,AL,AR,BD,BL,BR,b,Node);
	
	for(i=1;i<=NT;i++)
	{
		GS(AD,AL,AR,BD,BL,BR,x,xx,b,Node,eps,bc,D0,D1,N0,N1);
		for(j=0;j<Node;j++) x[j]=xx[j];
		
		if(i%1000==0)
		{
			cout<<i<<endl;
			for(j=0;j<Node;j++) fk<<x[j]<<" ";
			fk<<endl;
		}
	}
	
	return 0;
}


以下が計算結果の図です。初期に与えた三角形の熱の分布が時間が経過するにつれてどんどんなまっていく様子が観察できます。両端を冷やしているので棒がどんどん冷めていくのです。もっと時間進行させると領域全体で熱の分布が0に収束していくのが観察できます。


f:id:mutsumunemitsutan:20171121215655p:plain:w450

ゆくゆくはこのブログだけで有限要素法やいろいろな語学を学べる、self-contained(必要なものが全てそろった)もなのにしていきたいと考えています。なので、記事を読んでいて説明が不足しているなと思われたり、わかりづらいなと思われる箇所がありましたら報せて頂けると幸いです。

離散化の詳細はお待ちください。