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

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

【有限要素法】有限要素法のプログラミングを勉強するときにどの順序で学ぶのがよいか?【導入編】

有限要素法のプログラミングを勉強したいと思っている方はたくさんいると思います。私も最初はちんぷんかんぷんでした。Poisson方程式やら移流拡散方程式やらたくさん出てきます。有限要素法や数値計算の本を山ほど仕入れてきてうんうん唸っていました。そんな悪戦苦闘を通して、最近有限要素法がだんだんわかってきました。有限要素法を効率よく勉強するためには、以下の問題に対するプログラムを作成していくのが一番効率がよいと思います。


1次元Poisson方程式(1次元Laplace方程式でも可)

1次元定常移流拡散方程式

1次元非定常拡散方程式

1次元非定常移流拡散方程式

1次元非定常Burgers方程式

まずは1次元から始めるのがよいと思います。もちろん有限要素法が威力を発揮し始めるのは2次元以上ですがシンプルなものから勉強するべきです。

1次元Poisson方程式

そこで最初に1次元Poisson方程式を解くコードを作ります。これは拡散項のある定常問題を解けるようにするためです。ソース項のないLaplace方程式でも可です。ここで大規模な連立一次方程式の解法も学びます。私はGauss-Seidelを愛用していますが、Gaussの消去法も押さえておきましょう。TDMAでもokです。前者が反復法の代表で、後者が直接法の代表です。

1次元定常移流拡散方程式

この問題で拡散項と移流項が存在する定常問題の解き方を学びます。

1次元非定常拡散方程式

この問題では、拡散項だけが存在する非定常問題の解き方を学びます。陽解法、陰解法ともに作ります( \theta 法で離散化するとよいです)。

1次元非定常移流拡散方程式

次に非定常拡散方程式に移流項を足した問題の解き方を学びます。陽解法、陰解法ともに作ります。これが出来たら線型はもう大丈夫です。非線型問題を解くことは、線型問題を繰り返し解いていくのと同じなので線型問題が解けることは非常に重要です。

1次元非定常Burgers方程式

最後が1次元非定常Burgers方程式です。この問題では、非線型問題の解き方を学びます。Burgers方程式の場合は移流項が非線型項になります。将来Navier-Stokes方程式を解くときにも役立ちます。陽解法、陰解法ともに作ります。陰解法の場合は適当な反復法(Newton法とかPicard反復とか)が必要なのでそれも学びます。


導入なのでざっくりとしか書いていませんが流れはこんな感じです。このあとは風上化編に続きます(執筆中)。コードを書く際に参考になる文献の紹介や各問題の詳細な解説はまた別の記事でやります。

注 1次元移流方程式はリストに入れていませんが、これは風上化を導入した後に解くべきなのであとまわしになっています。

【読書リンク】読書猿とふろむだ

読書猿さんのブログ以来の大物です。これは本をしっかり読んでいる人の文章です。しかも並大抵の知識人ではないでしょう。読書猿さんのブログにより私の価値観は大いに影響を受けました。今度はふろむださんのブログが私の価値観に影響を与えるでしょう。

【数物リンク】熊野と学ぶ解析力学03【ラグランジュ方程式】

解析力学をやっているとラグランジアンというやつがいきなり出てきて面食らうのですが、「熊野と学ぶ解析力学03【ラグランジュ方程式】」は非常にうまくラグランジアンを導入しています。解析力学を勉強するのにとてもよくできた動画シリーズです。私も今見ています。

【Navier-Stokes方程式】MAC法によるNavier-Stokes方程式の離散化

今後キャビティー流れの計算をやりたいので、MAC法(Maker And Cell method)によるNavier-Stokes方程式の離散化について説明していきます。流速と圧力を直接未知数として計算する手法で陽解法です。

以下の内容は『流れ解析のための有限要素法』のpp.154-155を参考にしています。

流れ解析のための有限要素法入門

流れ解析のための有限要素法入門


簡単のため二次元直交座標系で考えます。まず保存型のNavier-Stokes方程式

 \displaystyle
\begin{eqnarray} 
\frac{\partial u}{\partial t} + \frac{\partial (u^2)}{\partial x} + \frac{\partial (uv)}{\partial y} &=& -\frac{\partial p}{\partial x} + \frac{1}{Re} \left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}\right) \\
\frac{\partial v}{\partial t} + \frac{\partial (uv)}{\partial x} + \frac{\partial (v)^2}{\partial y} &=& -\frac{\partial p}{\partial y} + \frac{1}{Re} \left( \frac{\partial^2 v}{\partial x^2} + \frac{\partial^2 v}{\partial y^2}\right)
\end{eqnarray}

です。連続式は

 \displaystyle \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0

です。Navier-Stokes方程式と連続式を以下のように時間に対して離散化します。

 \displaystyle
\begin{eqnarray} 
\frac{u^{n+1}-u^n}{\Delta t} + \frac{\partial (u^n)^2}{\partial x} + \frac{\partial (u^nv^n)}{\partial y} &=& -\frac{\partial p^n}{\partial x} + \frac{1}{Re} \left( \frac{\partial^2 u^n}{\partial x^2} + \frac{\partial^2 u^n}{\partial y^2}\right) \\
\frac{v^{n+1}-v^n}{\Delta t} + \frac{\partial (u^nv^n)}{\partial x} + \frac{\partial (v^n)^2}{\partial y} &=& -\frac{\partial p^n}{\partial y} + \frac{1}{Re} \left( \frac{\partial^2 v^n}{\partial x^2} + \frac{\partial^2 v^n}{\partial y^2}\right) \\
\frac{\partial u^{n+1}}{\partial x} + \frac{\partial v^{n+1}}{\partial y} &=& 0
\end{eqnarray}

次に1番目と2番目の式をn+1ステップでの流速について解いてそれらを3番目の式に代入すると

 \displaystyle
\begin{eqnarray} 
\frac{\partial^2 p^n}{\partial x^2} + \frac{\partial^2 p^n}{\partial y^2} &=& \frac{1}{\Delta t}  \left( \frac{\partial u^n}{\partial x} + \frac{\partial v^n}{\partial y} \right) -\frac{\partial^2 (u^n)^2}{\partial x^2} -\frac{\partial^2 u^n v^n}{\partial x \partial y} -\frac{\partial^2 (v^n)^2}{\partial y^2} \\
&+& \frac{1}{Re} \left[ \frac{\partial^2}{\partial x^2} \left( \frac{\partial u^n}{\partial x} + \frac{\partial v^n}{\partial y} \right) + \frac{\partial^2}{\partial y^2} \left( \frac{\partial u^n}{\partial x} + \frac{\partial v^n}{\partial y} \right)  \right]
\end{eqnarray}

を得ます。 \frac{\partial u^n}{\partial x} + \frac{\partial v^n}{\partial y} は連続式から0であるはずですが、数値計算上誤差が生じるので残してあります。しかし、 \frac{\partial u^n}{\partial x} + \frac{\partial v^n}{\partial y} はある程度小さいと考えられるのでその二階微分はほぼ0である

 \displaystyle \frac{\partial^2}{\partial x^2} \left( \frac{\partial u^n}{\partial x} + \frac{\partial v^n}{\partial y} \right) + \frac{\partial^2}{\partial y^2} \left( \frac{\partial u^n}{\partial x} + \frac{\partial v^n}{\partial y} \right) = 0

と考えることができます。よって

 \displaystyle
\frac{\partial^2 p^n}{\partial x^2} + \frac{\partial^2 p^n}{\partial y^2} = \frac{1}{\Delta t}  \left( \frac{\partial u^n}{\partial x} + \frac{\partial v^n}{\partial y} \right) -\frac{\partial^2 (u^n)^2}{\partial x^2} -\frac{\partial^2 u^n v^n}{\partial x \partial y} -\frac{\partial^2 (v^n)^2}{\partial y^2}

を得ます。これが圧力に関するポワソン方程式です。

計算手順は以下のように1と2の繰り返しとなります。

1. nステップでの流速  \boldsymbol{(u^n, v^n)} を用いて圧力のポワソン方程式を解きnステップでの圧力  \boldsymbol{p^{n+1}} を得る
2. nステップでの流速と圧力  \boldsymbol{(u^n, v^n, p^n)} を用いて離散化したNavier-Stokes方程式を解きn+1ステップでの流速  \boldsymbol{(u^{n+1}, v^{n+1})} を得る

空間方向の離散化には、保存型をそのまま用いれば有限体積法を、非保存型を用いれば有限要素法や差分法を用いることができます。

【有限体積法】1次元移流方程式を風上差分で解く C++コード付き

今回は1次元移流方程式を有限体積法(風上差分)で解いていきます。C++コード付きです。

1次元移流方程式とは

 \displaystyle \frac{\partial u}{\partial t} + V\frac{\partial u}{\partial x} =0, \quad x \in [0, L]

のような一階の偏微分方程式のことをいいます。ここで、 u は溶質の濃度、 V は風などによる移流速度です。

移流方程式はその名前の通り、初期の波形(プロファイル)が  V の速度で波形を変えずに移動していきます。移流方程式はこのように非常に直観的な解を持つのですが、これを数値計算で扱うのは非常に難しいです。いかに移流方程式をきとんと解くか、これがCFD(Computational Fluid Dynamics, 数値流体力学)における究極的な問いです。移流方程式をきちんと解くために今まで多くの労力がさかれてきました。なので、CFDを学ぼうと考えている方は移流方程式の数値解法をきちんと勉強することをおすすめします。今までにたくさんのスキームが移流方程式を解くために作られてきたのです。私も今後そうしたスキーム達を紹介していきます。

さて、風上差分というのは、特性曲線の向き(情報の向き)を考えることによって移流方程式を離散化しようという試みです。移流速度が正の場合は波形は左からくるので左のステンシルの値を用いるようにします。移流速度が負の場合はその逆です。詳しくは別記事で説明します。お待ちください。

境界条件は2パターンあります。すなわち、

Case 1,  \boldsymbol{x=0} でDirichlet条件を課す
 u(0)の値を指定)
Case 2,  \boldsymbol{x=L} でDirichlet条件を課す
 u(L) の値を指定)

です。それぞれコード上ではbc=1、bc=2に対応しています。特性曲線の向きを考えるとわかるのですが、 \boldsymbol{x=0} でDirichlet条件を課すのは移流速度が正の場合です。直観的には波が境界の外から内部へ入って来る様子を表しています。なので、 \boldsymbol{x=L} では境界条件は必要ありません。これが移流方程式の境界条件の独特な部分です。逆に移流速度が負の場合は x=L でDirichlet条件を課します。

今回の計算条件は、移流速度が正  V=0.1 で、分割数は400、時間刻みは0.001、領域の広さは  L=1、初期条件は領域の左端の矩形波(長方形の波)です。

以下では2つコードを示します。絶対値を使うものとC++の組み込み関数maxを使うコードです。どちらを使っても大丈夫です。絶対値を使う方法は離散化方程式を理論的に調べる際によく使われます。一方maxを使うとプログラムがとても簡単に書け、読みやすくなります。絶対値はこちらの記事を見てください。

まずは絶対値のコードです。

#include <iostream>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <string>
#include <sstream>

using namespace std;


inline void solve(double x[], int &Node, double &dt, double &dx, double &V)
{
	int i;
	
	for(i=1;i<Node-1;i++)
	{
		x[i]=x[i]-dt/dx*(-(V+abs(V))/2.0*x[i-1]+abs(V)*x[i]+(V-abs(V))/2.0*x[i+1]);
	}

}

inline void boundary(int &bc, double &D0, double x[], double &V, int &Node, double &dt, double &dx)
{
	if(bc==1)
	{
		x[0]=D0;
		x[Node-1]=x[Node-1]-dt/dx*V*(x[Node-1]-x[Node-2]);
	}
	if(bc==2)
	{
		x[0]=x[0]-dt/dx*V*(x[1]-x[0]);
		x[Node-1]=D0;
	}
}

inline void text(int &i,double x[], double &dx, int &Node)
{
	int j;
	
	stringstream ss;
	string name;
	ofstream fo;
	
	ss<<i;
	name=ss.str();
	name="Answer_" + name + ".txt";
	fo.open(name.c_str ());
	
	for(j=0;j<Node;j++)
	{
		fo<<dx*float(j)<<" "<<x[j]<<endl;
	}
}

int main()
{
	int i;
	int Partition=400;
	int Node=Partition+1;
	double LL=1.0;
	double dx=LL/(Node-1);
	double dt=0.001;
	double NT=7000;
	double eps=pow(2.0,-50);
	double V=0.1;
	int bc=1;//bc=1 left Dirichlet (V>=0), bc=2 right Dirichlet (V<0)
	double D0=0.0;//for Dirichlet
	double *x=new double[Node];
	
	ofstream fk;
	fk.open("Answer_0.txt");
	
	//initial condition//
	for(i=0;i<Node;i++)
	{
		x[i]=0.0;
		if((i>=40)&&(i<=80)) x[i]=1.0;
		fk<<dx*float(i)<<" "<<x[i]<<endl;
	}
	
	for(i=1;i<=NT;i++)
	{
		solve(x,Node,dt,dx,V);
		boundary(bc,D0,x,V,Node,dt,dx);
		
		if(i%1000==0)
		{
			cout<<i<<endl;
			text(i,x,dx,Node);
		}
	}
	
	delete[] x;
	
	return 0;
}

次にmaxを用いるコードです。

#include <iostream>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <string>
#include <sstream>

using namespace std;


inline void solve(double x[], int &Node, double &dt, double &dx, double &V)
{
	int i;
	
	for(i=1;i<Node-1;i++)
	{
		x[i]=x[i]-dt/dx*(-max(V,0.0)*x[i-1]+(max(V,0.0)+max(-V,0.0))*x[i]-max(-V,0.0)*x[i+1]);
	}

}

inline void boundary(int &bc, double &D0, double x[], double &V, int &Node, double &dt, double &dx)
{
	if(bc==1)
	{
		x[0]=D0;
		x[Node-1]=x[Node-1]-dt/dx*V*(x[Node-1]-x[Node-2]);
	}
	if(bc==2)
	{
		x[0]=x[0]-dt/dx*V*(x[1]-x[0]);
		x[Node-1]=D0;
	}
}

inline void text(int &i,double x[], double &dx, int &Node)
{
	int j;
	
	stringstream ss;
	string name;
	ofstream fo;
	
	ss<<i;
	name=ss.str();
	name="Answer_" + name + ".txt";
	fo.open(name.c_str ());
	
	for(j=0;j<Node;j++)
	{
		fo<<dx*float(j)<<" "<<x[j]<<endl;
	}
}

int main()
{
	int i;
	int Partition=400;
	int Node=Partition+1;
	double LL=1.0;
	double dx=LL/(Node-1);
	double dt=0.001;
	double NT=7000;
	double eps=pow(2.0,-50);
	double V=0.1;
	int bc=1;//bc=1 left Dirichlet (V>=0), bc=2 right Dirichlet (V<0)
	double D0=0.0;//for Dirichlet
	double *x=new double[Node];
	
	ofstream fk;
	fk.open("Answer_0.txt");
	
	//initial condition//
	for(i=0;i<Node;i++)
	{
		x[i]=0.0;
		if((i>=40)&&(i<=80)) x[i]=1.0;
		fk<<dx*float(i)<<" "<<x[i]<<endl;
	}
	
	for(i=1;i<=NT;i++)
	{
		solve(x,Node,dt,dx,V);
		boundary(bc,D0,x,V,Node,dt,dx);
		
		if(i%1000==0)
		{
			cout<<i<<endl;
			text(i,x,dx,Node);
		}
	}
	
	delete[] x;
	
	return 0;
}

計算結果を見てみましょう。

f:id:mutsumunemitsutan:20181111231051p:plain:w500

初期条件である左端の矩形波が時間がすすむにつれて右へと流されていきます。ただし、注意して頂きたいのは、厳密解だと初期波形がそのまま形を保ったまま流されていくはずですが、数値計算ではだんだんと角がとれてなまってきている点です。これは数値拡散と呼ばれています。風上差分は非常に安定して計算を行うことができるのですが、精度が1次しかなく、数値拡散が大きすぎます。そのため、この欠点を克服するために今までたくさんのスキームが作られてきました。

離散化の詳細は後でまとめます。

【有限体積法】1次元定常移流拡散方程式をハイブリッド法で解く C++コード付き

今回は1次元定常移流拡散方程式を有限体積法(ハイブリッド法)で解いていきます。C++コード付きです。

1次元定常移流拡散方程式とは

 \displaystyle V\frac{du}{dx} -D\frac{d^2u}{dx^2}=0, \quad x \in [0, L]

のような二階の常微分方程式のことをいいます。ここで、 u は溶質の濃度、 V は風などによる移流速度、 D は拡散の効果を表す拡散係数と解釈することができます。ある溶質が移流や拡散のもとでどのような分布をとるか求めるのがこの問題です。

前回は中心差分で解きましたが、移流が卓越する場合中心差分は不安定になるので対策が必要です。

今回はハイブリッド法を用います。ハイブリッド法とは簡単に言うと、移流が大きくない場合は中心差分を使い、移流が大きい場合は風上差分を使う方法のことです。境界条件は前回と同じく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法またはTDMAを用いています。好きなほうを選べます。詳しくは以下の記事を参考にして下さい。

コードです。

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

using namespace std;

//Gauss-Seidel//
inline void GS(double D[],double L[],double R[],double x[],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];
		
		for(i=0;i<Node;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]);
		}
		
		if(eps>max) break;
	}
}

//TDMA//a:diagonal,c:left,b:right,
inline void tdma(double a[], double c[], double b[], double d[], double x[], int size)
{
	int i;
	double *P=new double[size];
	double *Q=new double[size];
	
	//first step//
	P[0]=-b[0]/a[0];
	Q[0]=d[0]/a[0];
	
	//second step//
	for(i=1;i<size;i++)
	{
		P[i]=-b[i]/(a[i]+c[i]*P[i-1]);
		Q[i]=(d[i]-c[i]*Q[i-1])/(a[i]+c[i]*P[i-1]);
	}
	
	//third step, backward//
	x[size-1]=Q[size-1];

	for(i=size-2;i>-1;i=i-1)
	{
		x[i]=P[i]*x[i+1]+Q[i];
	}
	
	delete [] P,Q;
}

inline void mat(double D[],double L[],double R[],double b[],double f[], int &Node, double &dx, double &V, double &Diff)
{
	int i;
	for(i=1;i<Node;i++)
	{
		//Diffusion//
		L[i]=-max(V,max(Diff/dx+V/2.0,0.0));
		D[i]=max(V,max(Diff/dx+V/2.0,0.0))+max(-V,max(Diff/dx-V/2.0,0.0));
		R[i]=-max(-V,max(Diff/dx-V/2.0,0.0));
		
		//Source//
		b[i]=-f[i]*dx;
	}
}

inline void boundary(int &bc,double &D0,double &D1,double &N0,double &N1,double D[],double L[],double R[],double b[],int &Node, double &dx, double f[], double &V, double &Diff)
{
	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;
		D[0]=max(V,max(Diff/dx+V/2.0,0.0))-V;
		R[0]=-max(-V,max(Diff/dx-V/2.0,0.0));
		b[0]=-f[0]*dx/2.0-N0*Diff;
	}
	if(bc==3)
	{
		D[0]=1.0;R[0]=0.0;b[0]=D0;
		L[Node-1]=-max(V,max(Diff/dx+V/2.0,0.0));
		D[Node-1]=max(-V,max(Diff/dx-V/2.0,0.0))+V;
		b[Node-1]=-f[Node-1]*dx/2.0+N1*Diff;
	}
}

int main()
{
	int i;
	int Partition=100;
	int Node=Partition+1;//the number of cells is equal to the number of nodes
	double LL=1.0;
	double dx=LL/(Node-1);
	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=1.0; 
	double N0=5.0;
	double N1=3.0;
	double Diff=0.01;
	double V=0.1;
	double *b = new double[Node];
	double *x = new double[Node];
	double *f = new double[Node];
	double *D = new double[Node];
	double *L = new double[Node];
	double *R = new double[Node];

	for(i=0;i<Node;i++) 
	{
		b[i]=0.0;x[i]=0.0;f[i]=0.0;D[i]=0.0;L[i]=0.0;R[i]=0.0;
	}
	
	ofstream fk;
	fk.open("answer.txt");
	
	mat(D,L,R,b,f,Node,dx,V,Diff);
	boundary(bc,D0,D1,N0,N1,D,L,R,b,Node,dx,f,V,Diff);
	
	//Gauss-Seidel or TDMA//
	//GS(D,L,R,x,b,Node,eps);
	tdma(D,L,R,b,x,Node);
	
	for(i=0;i<Node;i++) cout<<x[i]<<endl;
	for(i=0;i<Node;i++) fk<<dx*i<<" "<<x[i]<<endl;
	
	delete[] b,x,f,D,L,R;

	return 0;
}

計算結果を見てみましょう。

f:id:mutsumunemitsutan:20181110130934p:plain:w500

中心差分の場合の結果と似ていますが、それはセル数が十分にあるからです。今後セル数が少ない場合の他の手法との比較をしたいと思います。

【有限体積法】1次元定常移流拡散方程式を風上差分で解く C++コード付き

今回は1次元定常移流拡散方程式を有限体積法(風上差分)で解いていきます。C++コード付きです。

1次元定常移流拡散方程式とは

 \displaystyle V\frac{du}{dx} -D\frac{d^2u}{dx^2}=0, \quad x \in [0, L]

のような二階の常微分方程式のことをいいます。ここで、 u は溶質の濃度、 V は風などによる移流速度、 D は拡散の効果を表す拡散係数と解釈することができます。ある溶質が移流や拡散のもとでどのような分布をとるか求めるのがこの問題です。

前回は中心差分で解きましたが、移流が卓越する場合中心差分は不安定になるので対策が必要です。

今回は風上差分を用います。境界条件は前回と同じく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法またはTDMAを用いています。好きなほうを選べます。詳しくは以下の記事を参考にして下さい。

以下では2つコードを示します。絶対値を使うものとC++の組み込み関数maxを使うコードです。どちらを使っても大丈夫です。絶対値を使う方法は離散化方程式を理論的に調べる際によく使われます。一方maxを使うとプログラムがとても簡単に書け、読みやすくなります。絶対値はこちらの記事を見てください。

まずは絶対値のコードです。

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

using namespace std;

//Gauss-Seidel//
inline void GS(double D[],double L[],double R[],double x[],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];
		
		for(i=0;i<Node;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]);
		}
		
		if(eps>max) break;
	}
}

//TDMA//a:diagonal,c:left,b:right,
inline void tdma(double a[], double c[], double b[], double d[], double x[], int size)
{
	int i;
	double *P=new double[size];
	double *Q=new double[size];
	
	//first step//
	P[0]=-b[0]/a[0];
	Q[0]=d[0]/a[0];
	
	//second step//
	for(i=1;i<size;i++)
	{
		P[i]=-b[i]/(a[i]+c[i]*P[i-1]);
		Q[i]=(d[i]-c[i]*Q[i-1])/(a[i]+c[i]*P[i-1]);
	}
	
	//third step, backward//
	x[size-1]=Q[size-1];

	for(i=size-2;i>-1;i=i-1)
	{
		x[i]=P[i]*x[i+1]+Q[i];
	}
	
	delete [] P,Q;
}

inline void mat(double D[],double L[],double R[],double b[],double f[], int &Node, double &dx, double &V, double &Diff)
{
	int i;
	for(i=1;i<Node;i++)
	{
		//Diffusion//
		L[i]=-Diff*1.0/dx;
		D[i]=Diff*2.0/dx;
		R[i]=-Diff*1.0/dx;
		
		//Advection//
		L[i]+=-(V+abs(V))/2.0;
		D[i]+=abs(V);
		R[i]+=(V-abs(V))/2.0;
		
		//Source//
		b[i]=-f[i]*dx;
	}
}

inline void boundary(int &bc,double &D0,double &D1,double &N0,double &N1,double D[],double L[],double R[],double b[],int &Node, double &dx, double f[], double &V, double &Diff)
{
	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;
		D[0]=(-V+abs(V))/2.0+Diff*1.0/dx;R[0]=(V-abs(V))/2.0-Diff*1.0/dx;b[0]=-f[0]*dx/2.0-N0*Diff;
	}
	if(bc==3)
	{
		D[0]=1.0;R[0]=0.0;b[0]=D0;
		L[Node-1]=-(V+abs(V))/2.0-Diff*1.0/dx;D[Node-1]=(V+abs(V))/2.0+Diff*1.0/dx;b[Node-1]=-f[Node-1]*dx/2.0+N1*Diff;
	}
}

int main()
{
	int i;
	int Partition=100;
	int Node=Partition+1;//the number of cells is equal to the number of nodes
	double LL=1.0;
	double dx=LL/(Node-1);
	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=1.0; 
	double N0=5.0;
	double N1=5.0;
	double Diff=0.01;
	double V=0.1;
	double *b = new double[Node];
	double *x = new double[Node];
	double *f = new double[Node];
	double *D = new double[Node];
	double *L = new double[Node];
	double *R = new double[Node];

	for(i=0;i<Node;i++) 
	{
		b[i]=0.0;x[i]=0.0;f[i]=0.0;D[i]=0.0;L[i]=0.0;R[i]=0.0;
	}
	
	ofstream fk;
	fk.open("answer.txt");
	
	mat(D,L,R,b,f,Node,dx,V,Diff);
	boundary(bc,D0,D1,N0,N1,D,L,R,b,Node,dx,f,V,Diff);
	
	//Gauss-Seidel or TDMA//
	//GS(D,L,R,x,b,Node,eps);
	tdma(D,L,R,b,x,Node);
	
	for(i=0;i<Node;i++) cout<<x[i]<<endl;
	for(i=0;i<Node;i++) fk<<dx*i<<" "<<x[i]<<endl;

	delete[] b,x,f,D,L,R;

	return 0;
}


次にmaxのコードです。

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

using namespace std;

//Gauss-Seidel//
inline void GS(double D[],double L[],double R[],double x[],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];
		
		for(i=0;i<Node;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]);
		}
		
		if(eps>max) break;
	}
}

//TDMA//a:diagonal,c:left,b:right,
inline void tdma(double a[], double c[], double b[], double d[], double x[], int size)
{
	int i;
	double *P=new double[size];
	double *Q=new double[size];
	
	//first step//
	P[0]=-b[0]/a[0];
	Q[0]=d[0]/a[0];
	
	//second step//
	for(i=1;i<size;i++)
	{
		P[i]=-b[i]/(a[i]+c[i]*P[i-1]);
		Q[i]=(d[i]-c[i]*Q[i-1])/(a[i]+c[i]*P[i-1]);
	}
	
	//third step, backward//
	x[size-1]=Q[size-1];

	for(i=size-2;i>-1;i=i-1)
	{
		x[i]=P[i]*x[i+1]+Q[i];
	}
	
	delete [] P,Q;
}

inline void mat(double D[],double L[],double R[],double b[],double f[], int &Node, double &dx, double &V, double &Diff)
{
	int i;
	for(i=1;i<Node;i++)
	{
		//Diffusion//
		L[i]=-Diff*1.0/dx;
		D[i]=Diff*2.0/dx;
		R[i]=-Diff*1.0/dx;
		
		//Advection//
		L[i]+=-max(V,0.0);
		D[i]+=max(V,0.0)+max(-V,0.0);
		R[i]+=-max(-V,0.0);
		
		//Source//
		b[i]=-f[i]*dx;
	}
}

inline void boundary(int &bc,double &D0,double &D1,double &N0,double &N1,double D[],double L[],double R[],double b[],int &Node, double &dx, double f[], double &V, double &Diff)
{
	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;
		D[0]=max(V,0.0)-V+Diff*1.0/dx;R[0]=-max(-V,0.0)-Diff*1.0/dx;b[0]=-f[0]*dx/2.0-N0*Diff;
	}
	if(bc==3)
	{
		D[0]=1.0;R[0]=0.0;b[0]=D0;
		L[Node-1]=-max(V,0.0)-Diff*1.0/dx;D[Node-1]=max(-V,0.0)+V+Diff*1.0/dx;b[Node-1]=-f[Node-1]*dx/2.0+N1*Diff;
	}
}

int main()
{
	int i;
	int Partition=100;
	int Node=Partition+1;//the number of cells is equal to the number of nodes
	double LL=1.0;
	double dx=LL/(Node-1);
	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=1.0; 
	double N0=5.0;
	double N1=5.0;
	double Diff=0.01;
	double V=0.1;
	double *b = new double[Node];
	double *x = new double[Node];
	double *f = new double[Node];
	double *D = new double[Node];
	double *L = new double[Node];
	double *R = new double[Node];

	for(i=0;i<Node;i++) 
	{
		b[i]=0.0;x[i]=0.0;f[i]=0.0;D[i]=0.0;L[i]=0.0;R[i]=0.0;
	}
	
	ofstream fk;
	fk.open("answer.txt");
	
	mat(D,L,R,b,f,Node,dx,V,Diff);
	boundary(bc,D0,D1,N0,N1,D,L,R,b,Node,dx,f,V,Diff);
	
	//Gauss-Seidel or TDMA//
	//GS(D,L,R,x,b,Node,eps);
	tdma(D,L,R,b,x,Node);
	
	for(i=0;i<Node;i++) cout<<x[i]<<endl;
	for(i=0;i<Node;i++) fk<<dx*i<<" "<<x[i]<<endl;
	
	delete[] b,x,f,D,L,R;

	return 0;
}

計算結果を見てみましょう。

f:id:mutsumunemitsutan:20181110090706p:plain:w500

中心差分の場合の結果と似ていますが、それはセル数が十分にあるからです。今後セル数が少ない場合の中心差分と風上差分の比較をしたいと思います。