【有限要素法】1次元Fisher-KPP方程式（フィッシャーKPP方程式）をGalerkin法で解く　C++コード付き

今回は1次元Fisher-KPP方程式（フィッシャーKPP方程式）を有限要素法（Galerkin法）で解いていきます。C++コード付きです。おそらくFisher-KPP方程式を有限要素法で解くコードが載っているネット上唯一のページです。Fisher-KPP方程式は半線型反応拡散方程式であり、生物が生息域を拡大している様子を表していると解釈することができます。

フィッシャーの方程式 - Wikipedia

Fisher's equation - Wikipedia

侵入・伝播と拡散方程式 (シリーズ・現象を解明する数学)

作者: 二宮広和,三村昌泰,竹内康博,森田善久
出版社/メーカー: 共立出版
発売日: 2014/07/09
メディア: 単行本
この商品を含むブログ (2件) を見る

1次元Fisher-KPP方程式とは

$\displaystyle \frac{\partial u}{\partial t}-D\frac{\partial^2 u}{\partial x^2}=ru \left(1-\frac{u}{K} \right), \quad x \in [0, L$ ]

のような偏微分方程式のことをいいます。ここで、 $u(t,x)$ は生物の生息密度、 $D$ は生物の広がり具合（拡散）を表す拡散係数、 $r$ が生物の成長率、 $K$ が環境容量（その生物が棲める最大密度）と解釈することができます。初期条件は $f(x)$ で与えられます。 $u$ は時間方向と空間方向ともに変化するので二変数関数 $u(t,x)$ となっています。右辺が二次の非線型項になっています。

Fisher-KPP方程式は拡散方程式

$\displaystyle \frac{\partial u}{\partial t}-D\frac{\partial^2 u}{\partial x^2}=0$

とロジスティック方程式（Logistic方程式）

$\displaystyle \frac{du}{dt}=ru \left(1-\frac{u}{K} \right)$

を組み合わせたものであると考えることができます。拡散によって生物の空間的な広がり（移動など）を、ロジスティック方程式によって（その場所における）生物の成長をあらわしているのです。似た話として線型反応移流拡散方程式が移流拡散方程式とマルサスモデルの組み合わせとして考えることができる、というのがあります。

プログラムより前に、まずは計算結果から見ていきましょう。2パターン見せます。パターン1は初期条件がsin波の生息密度の場合で、境界で密度が0のパターンです。パターン2は初期条件が至る所0で、片側から生物が移動してくるパターンです。詳しい計算条件については以下のように設定しています。まず、領域は $x\in[0,1$ ]、要素数は100、節点数は101で等分割（ $\Delta x=1/100$ ）、時間刻みは $\Delta t=0.001$ 、拡散係数は $D=0.001$ 、 $r=1$ 、 $K=1$ 、 $\theta=0.5$ です。パターン1については計算終了時刻は $t=8$ 、境界条件はbc=1（両側Dirichlet条件）で $u(t,0)=0$ と $u(t,1)=0$ 、初期条件はsin波で、式で書くと $f(x)=0.1\sin(\frac{\pi x}{L})$ です。パターン2については計算終了時刻は $t=21$ 、境界条件はbc=1（両側Dirichlet条件）で $u(t,0)=1$ と $u(t,1)=0$ 、初期条件は領域全体において0、式で書くと $f(x)=0$ です。

それではパターン1の計算結果です。

f:id:mutsumunemitsutan:20181026221857p:plain:w600

初期時刻に存在していた生物が増殖しながら領域全体に広がっていく様子が確認できます。生息密度は環境容量である $K=1$ へ向かって増えていきます。

次にパターン2の計算結果です。

f:id:mutsumunemitsutan:20181026222836p:plain:w600

初期時刻にはまったく領域において生物がいませんが、左側からどんどん領域へ侵入してきます。このとき右方向へと進む進行波解があらわれているのがわかります。要するに形を変えない解のことです。Fisher-KPP方程式には進行波解が存在することが広く知られています。この辺の話は『侵入・伝播と拡散方程式』や『反応拡散方程式』に詳しく書いてあります。生息密度は環境容量である $K=1$ へ向かって増えていきます。

侵入・伝播と拡散方程式 (シリーズ・現象を解明する数学)

作者: 二宮広和,三村昌泰,竹内康博,森田善久
出版社/メーカー: 共立出版
発売日: 2014/07/09
メディア: 単行本
この商品を含むブログ (2件) を見る

反応拡散方程式

作者: 柳田英二
出版社/メーカー: 東京大学出版会
発売日: 2015/06/26
メディア: 単行本
この商品を含むブログを見る

さて、中身の話に戻ります。非線型項の処理としてはPicard反復を使って線型化したFisher-KPP方程式を内部反復を入れて解いています。具体的には非線型項 $ru(1-\frac{u}{K}$ をソース項として処理しています。詳しくはPicard反復に関しては別記事で書きます。境界条件は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法またはTDMAを用いています。好きなほうを選べます。詳しくは以下の記事を参考にして下さい。

空間方向の離散化はGalerkin法、時間方向は $\theta$ 法で離散化しており、今回は $\theta=0.5$ のCrank-Nicolson法を用いています。ちなみに、 $\theta=0$ のとき前進Euler法（陽解法、explicit scheme）に、 $\theta=1$ のときに後退Euler法（完全陰解法、implicit scheme）になります。これらについても今後まとめる予定です。

コードは1次元Burgers方程式のものをもとに作りました。

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

using namespace std;

inline void GS(double AD[],double AL[],double AR[],double x[],double xx[],double b[], int &Node, double &eps)
{
	int i,k;
	double in,max;
	
	for(k=0;k<100000;k++)
	{
		max=0.0;
		for(i=0;i<Node;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]);
		}
		
		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 AD[],double AL[],double AR[],double BD[],double BL[],double BR[],double b[],double f[], int &Ele, double &dx, double &Diff, double &dt, double &theta, double x[], int &Node)
{
	int i;
	
	//initialization//
	for(i=0;i<Node;i++)
	{
		AD[i]=0.0;AL[i]=0.0;AR[i]=0.0;BD[i]=0.0;BL[i]=0.0;BR[i]=0.0;b[i]=0.0;
	}
	
	for(i=0;i<Ele;i++)
	{
		//A//
		//teporal//
		AD[i]+=dx*2.0/6.0/dt;
		AR[i]+=dx*1.0/6.0/dt;
		AL[i+1]+=dx*1.0/6.0/dt;
		AD[i+1]+=dx*2.0/6.0/dt;
		
		//diffusion//
		AD[i]+=theta*Diff/dx;
		AR[i]+=theta*-Diff/dx;
		AL[i+1]+=theta*-Diff/dx;
		AD[i+1]+=theta*Diff/dx;
		
		//B//
		//temporal//
		BD[i]+=dx*2.0/6.0/dt;
		BR[i]+=dx*1.0/6.0/dt;
		BL[i+1]+=dx*1.0/6.0/dt;
		BD[i+1]+=dx*2.0/6.0/dt;
		
		//diffusion//
		BD[i]+=-(1.0-theta)*Diff/dx;
		BR[i]+=-(1.0-theta)*-Diff/dx;
		BL[i+1]+=-(1.0-theta)*-Diff/dx;
		BD[i+1]+=-(1.0-theta)*Diff/dx;
		
		//Source//
		f[i]=x[i]*(1.0-x[i]);
		f[i+1]=x[i+1]*(1.0-x[i+1]);
		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 AD[],double AL[],double AR[],double BD[],double BL[],double BR[],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;
	}
}

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,j,k;
	int Ele=100;
	int Node=Ele+1;
	double LL=1.0;
	double dx=LL/Ele;
	double dt=0.001;
	double NT=8000;
	double eps=pow(2.0,-50);
	double eps2=pow(10.0,-10);
	double Diff=0.001;
	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;
	double max;
	double *b=new double[Node];
	double *x=new double[Node];
	double *xx=new double[Node];
	double *f=new double[Node];
	double *AD=new double[Node];
	double *AL=new double[Node];
	double *AR=new double[Node];
	double *BD=new double[Node];
	double *BL=new double[Node];
	double *BR=new double[Node];
	double *xr=new double[Node];
	
	ofstream fk;
	fk.open("Answer_0.txt");
	
	//initial condition//
	for(i=0;i<Node;i++)
	{
		x[i]=0.1*sin(4.0*atan(1.0)/LL*float(i)*dx);
		fk<<dx*float(i)<<" "<<x[i]<<endl;
	}
	
	
	for(i=1;i<=NT;i++)
	{ 
		//preserve x^n//
		for(j=0;j<Node;j++) xr[j]=x[j];
		
		//Picard iteration//
		for(j=0;j<100000;j++)
		{
			mat(AD,AL,AR,BD,BL,BR,b,f,Ele,dx,Diff,dt,theta,x,Node);
			boundary(bc,AD,AL,AR,BD,BL,BR,b,Node);
			
			//for b//
			if(bc==1)
			{
				b[0]+=D0;
				b[Node-1]+=D1;
			}
			if(bc==2)
			{
				b[Node-1]+=D1;
				b[0]+=BD[0]*xr[0]+BR[0]*xr[1]-N0*Diff;
			}
			if(bc==3)
			{
				b[0]+=D0;
				b[Node-1]+=BL[Node-1]*xr[Node-2]+BD[Node-1]*xr[Node-1]+N1*Diff;
			}
			
			for(k=1;k<Node-1;k++)
			{
				b[k]+=BL[k]*xr[k-1]+BD[k]*xr[k]+BR[k]*xr[k+1];
			}
	
			//GS or TDMA//
			//GS(AD,AL,AR,x,xx,b,Node,eps);
			tdma(AD,AL,AR,b,xx,Node);
			
			max=0.0;
			for(k=0;k<Node;k++)
			{
				if(max<abs(x[k]-xx[k])) max=abs(x[k]-xx[k]);
			}
			for(k=0;k<Node;k++) x[k]=xx[k];
			
			if(eps2>max) break;
		}
		
		if(i%1000==0)
		{
			cout<<i<<endl;
			text(i,x,dx,Node);
		}
	}
	
	delete[] b,x,xx,f,AD,AL,AR,BD,BL,BR,xr;
	
	return 0;
}

Burgers方程式以外の非線型偏微分方程式としてFisher-KPP方程式を取り上げてみました。

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