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

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

【数値計算】一次元非定常拡散方程式を有限要素法で解く 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(必要なものが全てそろった)もなのにしていきたいと考えています。なので、記事を読んでいて説明が不足しているなと思われたり、わかりづらいなと思われる箇所がありましたら報せて頂けると幸いです。

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

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

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

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

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

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

前回の一次元Poisson方程式に移流行列を加えるだけなので簡単です。今回も二次元配列を使うコードと一次元配列だけで処理するコードの両方を公開します。境界条件は前回と同じく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,double &D,double &V)
{
	int i;
	for(i=0;i<Ele;i++)
	{
		//Diffusion//
		A[i][i]+=D/dx;A[i][i+1]+=-D/dx;
		A[i+1][i]+=-D/dx;A[i+1][i+1]+=D/dx;
		
		//Advection//
		A[i][i]+=-V/2.0;A[i][i+1]+=V/2.0;
		A[i+1][i]+=-V/2.0;A[i+1][i+1]+=V/2.0;
		
		//Source//
		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 D=0.01;
	double V=0.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=0.0;
	double N1=2.0;
	vector<double> b(Node,0.);
	vector<double> x(Node,0.);
	vector<double> f(Node,0.);
	vector<vector<double> > A(Node,(vector<double>(Node,0.)));
	
	ofstream fk;
	fk.open("answer.txt");
	
	mat(A,b,f,Ele,dx,D,V);
	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;
}


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

#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, double &V, double &Diff)
{
	int i;
	for(i=0;i<Ele;i++)
	{
		//Diffusion//
		D[i]+=Diff/dx;R[i]+=-Diff/dx;
		L[i+1]+=-Diff/dx;D[i+1]+=Diff/dx;
		
		//Advection//
		D[i]+=-V/2.0;R[i]+=V/2.0;
		L[i+1]+=-V/2.0;D[i+1]+=V/2.0;
		
		//Source//
		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);
	double Diff=0.01;
	double V=0.1;
	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=0.0;
	double N1=2.0;
	vector<double> b(Node,0.);
	vector<double> x(Node,0.);
	vector<double> f(Node,0.);
	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,V,Diff);
	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;
}


こちらの離散化も詳しく書く予定です。Galerkin法では移流が卓越する場合をきちんと扱えないので、移流方向を考慮した離散化「Petrov-Galerkin法」を使う必要があります。それも後程。書きたいこと、書かねばならぬことが山のようにあります。ひとつずつ着実に消化していきたいです。

【数値計算】一次元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 です。

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

【数値計算】Gauss-Seidel法(ガウス・ザイデル法)のC++コード

現在「【数値計算】Gauss-Seidel法(ガウス・ザイデル法)の説明 実践編」を準備中ですが、それに先立ってGauss-Seidel法のC++コードを公開します。私はコーディングの専門家ではないのであまり綺麗ではありませんが悪しからず。皆様の勉強に役立てて頂ければ幸いです。解いている問題は

 \left( \begin{array}{ccc} 3 & 2 & 1 \\ 1 & 4 & 1 \\ 2 & 2 & 5 \end{array} \right) \left( \begin{array}{c} x_1 \\ x_2 \\ x_3 \end{array} \right) = \left( \begin{array}{c} 10 \\ 12 \\ 21 \end{array} \right)

で、答えは
 \left( \begin{array}{c} x_1 \\ x_2 \\ x_3 \end{array} \right) = \left( \begin{array}{c} 1 \\ 2 \\ 3 \end{array} \right)

となります。実際に動かして遊んでみて下さい。

#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 &size)
{
	double value=0.;
	int i;
	for(i=0;i<size;i++) value+=a[i]*b[i];
	return value;
}

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

int main()
{
	int size=3;
	int i;
	double eps=pow(2.0,-50);
	double in,max;
	vector<double> b(size,0.);
	vector<double> x(size,0.);
	vector<vector<double> > A(size,(vector<double>(size,0.)));
	
	A[0][0]=3;A[0][1]=2;A[0][2]=1;
	A[1][0]=1;A[1][1]=4;A[1][2]=1;
	A[2][0]=2;A[2][1]=2;A[2][2]=5;
	b[0]=10;b[1]=12;b[2]=21;
	
	GS(A,x,b,size,eps);

	for(i=0;i<size;i++) cout<<x[i]<<endl;

        return 0;
}

記事を書いていく順番について

ある事柄について書こうとすると、それを説明するための前提知識についても書かなくてはならなくなって、さらにその前提知識を説明するための前提知識についても書かなくてはならなくなって…、ときりがないのでこれからはどんどん書いていくことにします。もちろん、後で前提知識に関する記事を書いた場合はリンクを張る等していきます。あんまりしっかりやろうとして何も書けない状態よりは、多少がたがたしていてもどんどん書いていきます。

【数値計算】Gauss-Seidel法(ガウス・ザイデル法)の説明 理論編

Gauss-Seidel法とは

連立一次方程式を解くための基本的な反復法です。反復法とは一回の計算で解を求めるのではなく、ある決まった手順を何度も繰り返しながら真の解へと近づいていく方法です。数値計算では有限桁までしか表現できないので、得られる解は誤差を含んだものになります。ちなみに、何故連立一次方程式を解くのかというと、有限要素法などで方程式を離散化すると最終的に連立一次方程式に帰着されるからです。

反復手順

まず解こうとする連立一次方程式を

 Ax=b

としましょう。ここで、  A n \times n の行列、 x b n 次元のベクトルです。 A b が既知で、 x が未知です。例えば、 3 \times 3 だと

 \left( \begin{array}{ccc} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{array} \right) \left( \begin{array}{c} x_1 \\ x_2 \\ x_3 \end{array} \right) = \left( \begin{array}{c} b_1 \\ b_2 \\ b_3 \end{array} \right)

です。この行列  A

 A=D+L+U

と分解します。ここで、 D は対角行列(対角成分以外は成分が0の行列)、 L は下三角行列(対角線の下にだけ非零の成分が入っている行列)、 U は上三角行列(対角線の上にだけ非零の成分が入っている行列)です。いかめしく書いていますが、具体例を見ればすぐ分かると思います。例えば、

 A=\begin{pmatrix} 10 & 1 & 1 \\ 2 & 10 & 1 \\2 & 2 & 10 \\ \end{pmatrix}

としましょう。このとき、 A は、

 A=\begin{pmatrix} 10 & 1 & 1 \\ 2 & 10 & 1 \\ 2 & 2 & 10 \\ \end{pmatrix}=\begin{pmatrix} 10 & 0 & 0 \\ 0 & 10 & 0 \\ 0 & 0 & 10 \\ \end{pmatrix}+\begin{pmatrix} 0 & 0 & 0 \\ 2 & 0 & 0 \\ 2 & 2 & 0 \\ \end{pmatrix}+\begin{pmatrix} 0 & 1 & 1 \\ 0 & 0 & 1 \\ 0 & 0 & 0 \\ \end{pmatrix}

のように分解できるよね、と言っているだけです。さて、この分解を最初の連立一次方程式に代入すると

 (D+L+U)x=b

となります。次に、 Ux を右辺に移項して

 (D+L)x=-Ux+b

です。最後に、この両辺に  D+L逆行列  (D+L)^{-1} を左からかけて

 x=-(D+L)^{-1}Ux+(D+L)^{-1}b

を得ます。これをもとにGauss-Seidel法の反復式を作ります。まず、 k ステップ目での  x の近似解を  x^k としましょう。次のステップ( k+1 ステップ目)での  x の近似解を先ほどの式を使って

 x^{k+1}=-(D+L)^{-1}Ux^k+(D+L)^{-1}b

のように決めることにします。これがGauss-Seidel法の反復式です。適当な初期値  x^0 からスタートした数列  x^0, x^1,\cdots は、もし収束するなら  Ax=b を満たす  x へと収束します。なぜなら、数列  x^0, x^1,\cdots x^* に収束したとすると

 x^{*}=-(D+L)^{-1}Ux^*+(D+L)^{-1}b

が成り立ちます。後はこの式を逆に変形していけば

 Ax^*=b

となるので  x=x^* です。

実際の数値計算は有限桁なので、誤差が一定値以下になったら計算を切り上げます。誤差の決め方ですが、自分で決めた誤差の下限  \epsilon (小さい値)に対して

 \max_{1 \leq i \leq n} |x_i^{k+1} - x_i^{k}|\ < \epsilon

となったら反復を止めることにします。この式の意味は、 k+1 ステップ目の近似解と  k ステップ目の近似解の差を取ったとき、その差が最大になるものが、あらかじめ定めておいた誤差の下限を下回ったら反復を止めるということです。つまり、近似解の値があまり更新されなくなったら反復を止めよう、と言っているだけです。実は他にも止める基準はあるのですがそれはまた今度。

長くなってきたので一度切りますが、これだけではきっと分からないと思うので、具体的な行列に対する例は次回書きます。

注 Gauss-Seidel法の収束については、結構大変なので別の記事で触れます。

参考

英語で学ぶ数値解析

英語で学ぶ数値解析

楽しい反復法

楽しい反復法

【語学学習】どうやって新しい外国語を読めるようにするか?

私の外国語遍歴?

私は今までに色々な言語に手を出して参りました。大学での第二外学国語のアラビア語に始まり、ラテン語、ドイツ語、デンマーク語、オランダ語、古典ギリシャ語、フランス語、チェコ語、ロシア語、スペイン語、イタリア語。しかし、いろいろやりすぎて何も身に着かない状態が続いていました(それでも様々な言語に触れることが出来てよかったと思っていますが)。そこで、まずはひとつの言語をある程度までしっかり勉強してから他の言語に手を出すことにしました。そこで始めたのがドイツ語です。

文法を身につけましょう

ドイツ語を始めるにあたって、まずは薄い参考書を用いてその言語の全体をおおまかに把握することから始めました。あまり分厚い本に最初から取り組むと挫折しやすいので、最初は薄い参考書から始めましょう。この参考書をとにかく一度読み通します。活用や単語は、こんなのがあるんだ、ぐらいでどんどん進んでいきます。読み終えたらすぐに二週目を始めます。二回目は単語や活用をちゃんと覚えようと努力します。今ならAnkiを使って覚えます。これで基礎文法は大丈夫だと思います。今まで私が語学学習を挫折してきたのはこの二週目の復習が無かったから、だと思います。一回読んだくらいだと人間どんどん忘れてしまいます。欲を言えば三週すれば万全でしょう。

単語もちゃんと覚えましょう

実は、今まで私が語学学習を挫折してきた理由はもう一つあるんです。それは、基礎的な単語を1000語以上覚える努力を完全に怠っていたのです。入門書を一周してだいたいのその言語の特徴を把握するだけで満足していたのです。それではいけません!単語集等を使って覚えましょう。千野栄一先生も著書の中で述べられているように、「語学で必要なのは語彙と文法」なのです。さらに、「基礎語を1000覚えてしまえば初級は卒業でその言語は無に戻ることはない」ということなのです。だから新しい外国語をやる時は、薄い文法書1冊と基礎語1000語を覚えることができる環境を作り出さないといけません。これは例えば、よいテキストや単語集、外国語に取り組む時間、そして気力ないしやる気です。

私のドイツ語は、文法書をやった後、"Mastering German Vocabulary"や『ドイツ基本語5000辞典』を使ってちゃんと単語を覚えたので(後者は10月に終わる予定)、辞書をひきひきWikipediaを読めるようになりました。これで初級は卒業です。初級を卒業したら他の外国語に手を出しても大丈夫です。今、私もフランス語をはじめています。

薄い参考書が終わったら、同じようなレベルの参考書を何冊か読む、のもオススメです。文法の復習にもなるし、場合によってはある本には載っていない文法事項を学習できたりします(入門書は意外とレベルがまちまちでどこまで載っている範囲が全然違う)。さらにここで未知語を拾っておくと1000語の壁を超えるのがかなり楽になります。

外国語上達法 (岩波新書 黄版 329)

外国語上達法 (岩波新書 黄版 329)

Mastering German Vocabulary: A Thematic Approach (Mastering Vocabulary)

Mastering German Vocabulary: A Thematic Approach (Mastering Vocabulary)

  • 作者: Veronika Schnorr,Martin Crellin,Adelheid Schnorr-Dummler,Gabriele Forst,Raymond Sudmeyer
  • 出版社/メーカー: Barrons Educational Series Inc
  • 発売日: 1995/08
  • メディア: ペーパーバック
  • クリック: 4回
  • この商品を含むブログ (3件) を見る
ドイツ基本語5000辞典

ドイツ基本語5000辞典

そして終わりのない中・上級へ

難しいのが中級からです。これはもうどんどん原文を読んでいくしかありません。なるべく自分の興味がある内容がよいです。Wikipediaはこのためにとても役に立ちます。あらかじめ日本語や英語でその記事を読んでおけば、ドイツ語の記事の内容が想像できます。類推もかなり効くようになります。しかも、専門性の高い記事ほど英語を訳しただけだったりするのでかなり読めます。さらに関連事項も読んでいくとよいです。Wikipediaを読みつつ、未知語をAnkiに放り込み、わからない文法事項は詳しめの文法書で調べて解決しましょう。やはり、読んでいて自然と出会う語彙がその人にとって本当に必要な語彙である、と思います。あと多少メジャーな言語だと中級向けの解釈の参考書や読本があるのでそれで学習しましょう。基本的に英語以外の言語の参考書はいつ絶版になってAmazonで高額で購入しないといけなくなるかわからないので積極的に入手していきましょう。私も将来やるであろう外国語の参考書を普段から集めています。趣味みたいなもんです。

注 新しい外国語は同時並行で2つ以上勉強しないほうがいいです。経験上どちらもポシャります。

参考