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

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

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

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

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

Fisher's equation - Wikipedia

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

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

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 へ向かって増えていきます。

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

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

反応拡散方程式

反応拡散方程式


さて、中身の話に戻ります。非線型項の処理としては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方程式を取り上げてみました。

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

【有限要素法】1次元Poisson方程式(ポワソン方程式)の解析解

1次元Poisson方程式

 \displaystyle \frac{d^2u}{dx^2}+f(x)=0, \quad x \in [0, L]

の解析解を求めます。 f=1 とし、境界条件 u(0)=0 \left. \frac{du}{dx} \right|_{x=L} = 0 とします。まず  f を右辺に移項して一回積分すると

 \displaystyle \frac{du}{dx}=-\int 1 dx + C_1 = -x + C_1

を得ます。ここで  C_1積分定数です。 \left. \frac{du}{dx} \right|_{x=L} = 0 より

 \displaystyle 0 = -L + C_1

です。すなわち  \displaystyle L = C_1 を得ます。これを代入して

 \displaystyle \frac{du}{dx}= -x + L

となります。次にもう一度積分すると

 \displaystyle u= - \int (x+L) dx +C_2 = -\frac{x^2}{2} + Lx + C_2

となります。 u(0)=0 より

 \displaystyle 0= C_2

です。これを代入して結局解は

 \displaystyle u= -\frac{x^2}{2} + Lx

です。


この解を用いて1次元Poisson方程式に対する有限要素法(Galerkin法)の精度を調べます。

【有限要素法】1次元Burgers方程式(バーガース方程式)をGalerkin法で解く C++コード付き

今回は1次元Burgers方程式(バーガース方程式)を有限要素法(Galerkin法)で解いていきます。C++コード付きです。1次元Burgers方程式は、前に書いた記事「有限要素法のプログラミングを勉強するときにどの順序で学ぶのがよいか?」で紹介した、有限要素法を学ぶ場合に五番目に解くべき方程式です。この問題を解くことで非線型の場合の扱いを学びます。ただし今回は風上化は入れていません。

1次元Burgers方程式とは

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

のような偏微分方程式のことをいいます。ここで、 u(t,x) は熱や溶質の濃度、 D は拡散の効果を表す拡散係数と解釈することができます。初期条件は  f(x) で与えられます。1次元非定常移流拡散方程式との違いは移流項の前にかかっている係数です。1次元非定常移流拡散方程式ではこの係数は定数でしたが、1次元Burgers方程式ではこの係数は解そのものとなっており、ここが非線型項になっています。 u は時間方向と空間方向ともに変化するので二変数関数  u(t,x) となっています。非線型項の処理としてはPicard反復(ピカード反復ないしはピカール反復、前者は英語読みで後者はフランス語読み。Picardはフランス人。)を使って線型化したBurgers方程式を内部反復を入れて解いています。ここで書くと長くなるので別記事で書きます。境界条件は前回と同じく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)になります。これらについても今後まとめる予定です。

計算条件については以下のように設定しています。まず、Burgers方程式を解く領域は  x\in[0,1]、要素数は100、節点数は101で等分割( \Delta x=1/100)、時間刻みは  \Delta t=0.001、計算終了時刻は  t=1、拡散係数は  D=0.01 \theta=0.5境界条件はbc=1(両側Dirichlet条件)で  u(t,0)=0 u(t,1)=0、初期条件はsin波で、式で書くと  f(x)=\sin(\frac{\pi x}{L}) です。初期時刻にsin波の分布を与え、両側での値を0に固定した場合に、波がある時間にどのような分布をとるか、という問題です。答えは"Answer_****.txt"(****にはステップ数が入る)に出力されます。100ステップごとに値が出力されます。一行目が  x 座標、二行目がその時刻における  u の分布となっています。

コードです。

#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;
	double VE;
	
	//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;
	}
	
	for(i=0;i<Ele;i++)
	{
		
		
		VE=(x[i]+x[i+1])/2.0;
		
		//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;
		
		//advection//
		AD[i]+=theta*-VE/2.0;
		AR[i]+=theta*VE/2.0;
		AL[i+1]+=theta*-VE/2.0;
		AD[i+1]+=theta*VE/2.0;
		
		//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;
		
		//advection//
		BD[i]+=-(1.0-theta)*-VE/2.0;
		BR[i]+=-(1.0-theta)*VE/2.0;
		BL[i+1]+=-(1.0-theta)*-VE/2.0;
		BD[i+1]+=-(1.0-theta)*VE/2.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=1000;
	double eps=pow(2.0,-50);
	double eps2=pow(10.0,-10);
	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;
	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]=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%100==0)
		{
			cout<<i<<endl;
			text(i,x,dx,Node);
		}
	}
	
	delete[] b,x,xx,f,AD,AL,AR,BD,BL,BR,xr;
	
	return 0;
}


以下が計算結果の図です。初期に与えたsin波の分布が時間が経過するにつれてどんどんなまりつつ右へと流されていく様子が観察できます。その際値が大きいところほど右へ動く速度が速く、だんだん波が切り立つようになっていく様子が観察できます。

f:id:mutsumunemitsutan:20181023204834p:plain

非線型」という用語の意味は「線型でない」なので非線形問題を一括りにするのは不可能です。ケースバイケースで対応しなければなりません。非線型偏微分方程式の一例としてFisher-KPP方程式があります。この場合はソース項が二次で非線型となっています。


1次元のBurgers方程式はCole-Hopf変換によって解析的に解ける場合があります。


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

【有限要素法】1次元非定常移流拡散方程式をGalerkin法で解く C++コード付き

今回は1次元非定常移流拡散方程式を有限要素法(Galerkin法)で解いていきます。C++コード付きです。1次元非定常移流拡散方程式は、前に書いた記事「有限要素法のプログラミングを勉強するときにどの順序で学ぶのがよいか?」で紹介した、有限要素法を学ぶ場合に四番目に解くべき方程式です。この問題は線型の場合の総まとめとなっています。1次元非定常移流拡散方程式が解けるようになれば基本的に何でも有限要素法で解けるようになります。なのでしっかり勉強してください。ただし今回は風上化は入れていません。

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

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

のような偏微分方程式のことをいいます。ここで、 u(t,x) は熱や溶質の濃度、 V は流速、 D は拡散の効果を表す拡散係数と解釈することができます。初期条件は  f(x) で与えられます。熱またはある溶質が時間が経過するにつれてどのような分布をとるか(移流・拡散していくか)求めるのがこの問題です。 u は時間方向と空間方向ともに変化するので二変数関数  u(t,x) となっています。前に説明した1次元移流拡散方程式を時間発展させたものにあたります。境界条件は前回と同じく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)になります。これらについても今後まとめる予定です。

計算条件については以下のように設定しています。まず、非定常移流拡散方程式を解く領域は  x\in[0,1]、要素数は100、節点数は101で等分割( \Delta x=1/100)、時間刻みは  \Delta t=0.001、計算終了時刻は  t=10、流速は  V=0.1、拡散係数は  D=0.01 \theta=0.5境界条件はbc=1(両側Dirichlet条件)で  u(t,0)=0 u(t,1)=0、初期条件は頂点の値が0.5の三角形(下の図参照)です。式で書くと  f(x)=\min(x,1-x) です。初期時刻に三角形のような濃度の分布を与え、棒の両側での濃度を0に固定し、右向きに流れがある場合に、濃度がある時間にどのような分布をとるか、という問題です。答えは"Answer_****.txt"(****にはステップ数が入る)に出力されます。1000ステップごとに値が出力されます。一行目が  x 座標、二行目がその時刻における  u の分布となっています。

コードです。

#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 &V, 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;
	}
	
	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;
		
		//advection//
		AD[i]+=theta*-V/2.0;
		AR[i]+=theta*V/2.0;
		AL[i+1]+=theta*-V/2.0;
		AD[i+1]+=theta*V/2.0;
		
		//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;
		
		//advection//
		BD[i]+=-(1.0-theta)*-V/2.0;
		BR[i]+=-(1.0-theta)*V/2.0;
		BL[i+1]+=-(1.0-theta)*-V/2.0;
		BD[i+1]+=-(1.0-theta)*V/2.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;
	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 V=0.1;
	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 *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];
	
	ofstream fk;
	fk.open("Answer_0.txt");
	
	//initial condition//
	for(i=0;i<Node;i++)
	{
		x[i]=min(dx*float(i),1.0-dx*float(i));
		fk<<dx*float(i)<<" "<<x[i]<<endl;
	}
	
	mat(AD,AL,AR,BD,BL,BR,b,f,Ele,dx,Diff,dt,theta,V,Node);
	boundary(bc,AD,AL,AR,BD,BL,BR,b,Node);
	
	for(i=1;i<=NT;i++)
	{
		//for b//
		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*Diff;
		}
		if(bc==3)
		{
			b[0]=D0;
			b[Node-1]=BL[Node-1]*x[Node-2]+BD[Node-1]*x[Node-1]+N1*Diff;
		}
		
		for(j=1;j<Node-1;j++)
		{
			b[j]=BL[j]*x[j-1]+BD[j]*x[j]+BR[j]*x[j+1];
		}

		//GS or TDMA//
		GS(AD,AL,AR,x,xx,b,Node,eps);
		//tdma(AD,AL,AR,b,xx,Node);
		for(j=0;j<Node;j++) x[j]=xx[j];
		
		if(i%1000==0)
		{
			cout<<i<<endl;
			text(i,x,dx,Node);
		}
	}
	
	delete[] b,x,xx,f,AD,AL,AR,BD,BL,BR;
	
	return 0;
}


以下が計算結果の図です。初期に与えた三角形の濃度の分布が時間が経過するにつれてどんどんなまりつつ右へと流されていく様子が観察できます。

f:id:mutsumunemitsutan:20181022213314p:plain:w500

どうですか、だんだん楽しくなってきたでしょう。私は楽しいです。

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

【数物リンク】2変数の微分積分学の基本定理(ライプニッツ則)

2変数の微分積分学の基本定理の証明が載っているリンクです。別名ライプニッツ則ともいいます。積分する領域が時間とともに変わる場合を扱うときに使います。例えばNavier-Stokes方程式から浅水流方程式を導く際に使います。よく忘れてしまうのでメモです。

https://www.chart.co.jp/subject/sugaku/suken_tsushin/38/38-6.pdf

【有限要素法】無限領域における拡散方程式の解析解

今回は無限領域における拡散方程式の解析解をフーリエ変換を使って導きたいと思います。松下泰雄『フーリエ解析』のpp.126-127を参照します。

フーリエ解析―基礎と応用

フーリエ解析―基礎と応用

解く拡散方程式は

 \displaystyle \frac{\partial u}{\partial \tau} = D\frac{\partial^2 u}{\partial y^2} in  (-\infty, \infty)

です。初期条件はある関数  f(y) で与えられるとします。無限領域なので境界条件は必要なく、初期条件だけです。

変数  y に対するフーリエ変換を考えます。


 
\displaystyle 
\begin{eqnarray} 

\mathcal{F} [u(\tau,y)] &=& U(\omega,\tau) \\
\mathcal{F} [f(y)] &=& F(\omega) 

\end{eqnarray}

拡散方程式の両辺をフーリエ変換すると

 \displaystyle \mathcal{F} \left[\frac{\partial u}{\partial \tau}\right] = \mathcal{F} \left[D\frac{\partial^2 u}{\partial y^2}\right]

となります。まず左辺から考えましょう。 y に対するフーリエ変換を考えているので時間に対する偏微分を以下のように外に出せます。

 \displaystyle \mathcal{F} \left[\frac{\partial u}{\partial \tau}\right] = \frac{\partial}{\partial \tau} \mathcal{F} [u] = \frac{d}{d \tau} U(\omega,\tau)

最後のイコールで偏微分が常微分に変わるのは、 U \tau にだけ依存すると考えているからです。次に右辺を考えましょう。フーリエ変換の公式

 \displaystyle \mathcal{F} \left[ f^{(n)}(x) \right] = (i\omega)^n F(\omega)

を使います(p.71)。この公式はn階微分した関数のフーリエ変換はもとの関数のフーリエ変換 (i\omega) をかけたものに等しい、と言っています。よって

 \displaystyle \mathcal{F} \left[D\frac{\partial^2 u}{\partial y^2}\right] = -D\omega^2 \mathcal{F} [u] = -D\omega^2 U(\omega,\tau)

となるので結局  \tau に対する常微分方程式

 \displaystyle \frac{d}{d \tau} U(\omega,\tau) = -D\omega^2 U(\omega,\tau)

を得ます。これは変数分離で簡単に解けて(線型一階常微分方程式です)

 \displaystyle U(\omega,\tau) = A(\omega) e^{-D\omega^2 \tau}

となります。ここで  A \tau に対する積分定数です。これは初期条件  F(\omega) を使って消します。 t=0 のときを考えると

 \displaystyle U(\omega, 0) = A(\omega) = F(\omega)

となるので

 \displaystyle U(\omega,\tau) = F(\omega) e^{-D\omega^2 \tau}

が解です。

次はこれをフーリエ逆変換してもとの世界に戻します。まず、ガウス関数  e^{-ax^2} a>0フーリエ変換を考えます(p.66)。

 \displaystyle \mathcal{F} \left[ e^{-ay^2} \right] = \sqrt{\frac{\pi}{a}} e^{-\frac{\omega^2}{4a}}

ここで  a=\frac{1}{4D\tau} とおくと上式は

 \displaystyle \mathcal{F} \left[ e^{-\frac{y^2}{4D\tau}} \right] = \sqrt{4D\tau\pi} e^{-D\omega^2 \tau}

となります。これを使うと  U の式は

 \displaystyle U(\omega,\tau) = \mathcal{F} [f(y)] \mathcal{F} \left[ \frac{1}{\sqrt{4D\tau\pi}} e^{-\frac{y^2}{4D\tau}} \right]

と2つのフーリエ変換の積でかけます。ここで畳み込みの公式(p.77)

 \displaystyle \mathcal{F} [f \ast g ] = \mathcal{F} [f] \mathcal{F} [g]

を使います。ただし

 \displaystyle \mathcal{F} [f \ast g (y)] = \int_{-\infty}^{\infty} f(s)g(x-s) ds

です。これを使うと  U の式は

 \displaystyle U(\omega,\tau) = \mathcal{F} \left[f(y) \ast \frac{1}{\sqrt{4D\tau\pi}} e^{-\frac{y^2}{4D\tau}} \right]

とかけます。この式をフーリエ逆変換すると

 \displaystyle u(\tau,y) = f(y) \ast \frac{1}{\sqrt{4D\tau\pi}} e^{-\frac{y^2}{4D\tau}}

です。積分でかくと

 \displaystyle u(\tau,y) = \frac{1}{\sqrt{4D\tau\pi}} \int_{-\infty}^{\infty} f(s)  e^{-\frac{(y-s)^2}{4D\tau}} ds

となり拡散方程式の解を得られました。解析解と初めに言いましたが、この積分を解析的に求められる場合以外は解析解にはならないのでご注意下さい。


ちなみにこの解は1次元非定常移流拡散方程式や1次元Burgers方程式の解析解を導く際に用いられます。

【クラシック】超有望ヴァイオリニスト、Chloe Chua

今回紹介するのはシンガポールの超有望ヴァイオリニスト、Chloe Chuaです。なんと12歳です。2018年のユーディ・メニューイン国際コンクールで優勝しました。とにかく動画を見て下さい。おったまげます。


ユーディ・メニューイン国際コンクールでのガラコンサート(記念演奏会)、曲目は「ヴィヴァルディの四季」より冬全楽章


ヴィエニャフスキの創作主題による華麗なる変奏曲」(Wieniawski, Variations on an Original Theme, Op. 15)


ユーディ・メニューイン国際コンクール、セミファイナル、曲目は30秒から「ベートーベンのヴァイオリンソナタ8番」より第1楽章(Beethoven, Violin Sonata No. 8 in G major, Op. 30, No. 3, 1st movement )、5分30秒から「ピアソラのタンゴの歴史」よりカフェ 1930(Piazzolla, Histoire du Tango, Café 1930)、13分45秒から「フバイのカルメンによる華麗な幻想曲」(Hubay, Carmen Fantasie Brillante)


堂々としたソロ、無伴奏で顕著にあらわれている高度なフレーズの解釈力、盤石な音程、まったく違うスタイルの曲を表現する表現力、とても12歳とは思えません。Youtubeで適当に自動再生でクラシックを流していてたまたま流れたのがChloeの演奏でした。最初はこれはかなり有望なほぼプロと言ってよいレベルの20代前半の演奏家だと思ったのですが、動画を見て12歳でここまで到達できるものなのかと愕然としました。才能と努力の両方があれどこれは難しいのではないでしょうか。普通にクラシックとして鑑賞に堪えるレベルなのです。むしろ好んで聞きたくなるほどです。今まで幾多の「神童」がいましたが、彼らは年齢の割に上手であるとか、難曲を一応弾けるといったレベルでした。しかし、Chloeは違います。彼女はプロの演奏家と比べて遜色無いと私は思います。きっと彼女は次の時代の巨匠になります。楽しみです。

何かで読みましたが、クラシックの本場であるヨーロッパではクラシックの人気は下がりつつあり人も集まらなくなっているが、その代わりアジアなどの新興国から人が参入してきているそうです。これは日本で相撲をやる人が減りモンゴルから有望な力士がやってくる構図と同じだそうです。


追記
最新の動画があがってきました。


サラサーテカルメン幻想曲 」(Sarasate, Carmen Fantasy, Op. 25)

芸に磨きがかかっています!まったく危なげなく弾きこなしていますね。高音、フラジオレット、重音が気持ちよいです。カルメン幻想曲には多くの弱く弾く箇所が出てきますが、メリハリをつけて弱い箇所を弱く弾けるのが素晴らしいです。サラサーテカルメン幻想曲を聞くとき、私はいつも重音が連続する箇所に注目して聞いています。上記の動画だと10:09からです。chloeはきちんと弾きこなしています。プロを自称する方でもここが弾けていない例がたくさんあります。音程が合っていなかったり音がつながっていたり。Youtubeで適当に調べてみてください。うんざりします。彼女の演奏はギル・シャハムの演奏を思わせる快演です!!!いずれギル・シャハムに匹敵するヴァイオリニストになるでしょう。



メンデルスゾーンのヴァイオリン協奏曲」(Mendelsohn, Violin Conserto in E minor, Op. 64)

安定感と丁寧さはヒラリーハーンを思わせます。


www.youtube.com
ミルシテインパガニーニアーナより」(Milstein, Paganiniana)

パガニーニアーナは名ヴァイオリニストであるミルシテインパガニーニの24のカプリースをもとに作曲した無伴奏の超難曲です。今回は一部だけです。是非とも全部聴きたいですね。


参考

サラサーテカルメン幻想曲 」(Sarasate, Carmen Fantasy, Op. 25)
ギル・シャハムによる演奏、疾走感が凄まじい。例の箇所は8:52より。この速さで完璧に弾いている。サラサーテカルメン幻想曲の一つの完成形であると言いたい。欲を言えばこれを超える演奏を聴いてみたい。


ミルシテインパガニーニアーナ」(Milstein, Paganiniana)
本人による演奏。何度聞いても聴き飽きない。一生聴ける。難しいパッセージなのに主旋律がきちんと聴こえる。すごい。