差分スキーム物理化学生物現象には微分方程式でモデル化される例が多いモデルを使って現実の現象をコンピュータ上で再現することをシミュレーション ( 数値シミュレーションコンピュータシミュレーション ) と呼ぶそのためには微分方程式をコンピュータ上で計算できる数値スキームで近似することが必要

Size: px

Start display at page:

Download "差分スキーム物理化学生物現象には微分方程式でモデル化される例が多いモデルを使って現実の現象をコンピュータ上で再現することをシミュレーション ( 数値シミュレーションコンピュータシミュレーション ) と呼ぶそのためには微分方程式をコンピュータ上で計算できる数値スキームで近似することが必要"

えみこいまる
5 years ago
Views:

1 差分スキーム物理化学生物現象には微分方程式でモデル化される例が多いモデルを使って現実の現象をコンピュータ上で再現することをシミュレーション ( 数値シミュレーションコンピュータシミュレーション ) と呼ぶそのためには微分方程式をコンピュータ上で計算できる数値スキームで近似することが必要になるその一つの方法が微分方程式を差分方程式におき直すことである微分方程式の差分化次の 1 次元境界値問題を考える格子上の関数区間 (0,1) を分割して各格子点上での関数値を設定する区間 (0,1) をn 分割すれば格子点数はn+1 で格子間隔 h= 格子点の位置格子点での u の値 u()= 格子点での関数値 1+2 を定義できるまた次の値を定義する u( u 1 次微分 u'(x), 2 次微分 u''(x) の差分を導出する Taylor 展開による方法 1 次微分 u'(x) の差分 Maple のtaylor コマンドを使ってu(x) をTaylor 展開する ua:=taylor(u(x+h),h,2); ub:=taylor(u(x-h),h,2); 前進差分 u(x+h) のTaylor 展開より (1) (2) 後退差分 u(x-h) のTaylor 展開より

2 A uaa:=taylor(u(x+h),h,3); ubb:=taylor(u(x-h),h,3); (3) (4) 中心差分 uaa-ubb を計算して u'(x) 2 次微分 u''(x) の差分 uaa+ubb を計算して 2 次微分 u''(x) の差分を求めることができる差分スキームの構成例 11.1 の差分スキームを求める i=1, $$$n-1 これを行列で表現すると (n-1) (n-1) 行列として与えられた n に対して連立 1 次方程式 Au=b を解くことによって境界値問題の近似解として差分解を求めることができる

3 差分スキームのプログラミングステップ 1 格子点の生成例えば区間 (0,1) をn=6 分割する n:=7; gridpoints:=n; h:=evalf(1/gridpoints,7); x:=vector(1..n+1): for i from 1 to n+1 do: x[i]:=h*(i-1); end do; (5) (6) (7) (8) ステップ 2 係数行列の設定 A:=Matrix(1..n-1,1..n-1): for i from 1 to n-2 do: A[i,i+1]:=-1.0/(h*h): A[i+1,i]:=-1.0/(h*h): for i from 1 to n-1 do: A[i,i]:=2.0/(h*h)+x[i+1]: A; (9) ステップ 3 ベクトル b の設定 b:=vector(1..n-1);

4 (10) for i from 1 to n-1 do: b[i]:=( *x[i+1]-x[i+1]^2)*e end do; (11) ステップ 4 境界条件の設定 b[1]:=b[1]+1.0/(h*h); (12) ステップ 5 行列の解法 with(linearalgebra): u:=linearsolve(a,b); (13) ステップ 6 境界値の付加 uu:=vector(1..n+1): for i from 1 to n-1 do uu[i+1]:=u[i]; uu[1]:=1.0: uu[n+1]:=0.0: f1:=plot(x,uu):

5 f2:=plot((1-t)*exp(t),t= ,colo with (plots): display(f1,f2); t 非定常な問題 : 拡散現象時間発展する非定常な問題として拡散現象を考える 1 次元初期値問題 : 拡散方程式 ) 格子上の関数空間に対して区間 (-K,K) を分割し時間に対して区間 (0,T) を分割して各格子点上での関数値を設定する区間 (-K,K) をn 分割し区間 (0,T) をm 分割する格子点数はn+1 で格子間隔 h= 格子点の位置

6 A を定義できる時刻格子点での u の値 u(,)= 格子点での関数値時間の差分スキーム陽的スキーム u = 陰的スキーム u = とすると陽的スキーム陰的スキーム差分スキームのプログラミング ( 陰的スキームと陽的スキーム ) ステップ 1 格子点の生成例えば空間の区間 (-10,10) をn=50 分割, 時間の区間 (0,1.6) を10 分割する D 1とする n:=50; gridpoints:=n; (14) (15)

7 h:=evalf(20/gridpoints,7); deltat:=0.16; x:=vector(1..n-1): for i from 1 to n-1 do: x[i]:=-10.0+h*i; (15) (16) (17) ステップ 2 係数行列の設定 E:=Matrix(n-1,n-1,shape=identity): A:=Matrix(1..n-1,1..n-1): for i from 1 to n-2 do: A[i,i+1]:=-1.0/(h*h): A[i+1,i]:=-1.0/(h*h): for i from 1 to n-1 do: A[i,i]:=2.0/(h*h): B:=E+deltat*A: C:=E-deltat*A: A: ステップ 3 初期値の設定 u0:=vector(1..n-1): uk:=vector(1..n-1): vk:=vector(1..n-1): uu:=vector(1..n-1): for i from 1 to n-1 do: u0[i]:=exp(-x[i]*x[i]): uk:=u0: vk:=u0: with(linearalgebra): m:=10; k:=1; ステップ 4 解法 (18) (19) for k from 1 to m do: uu:=uk: v:=c.vk: u:=linearsolve(b,uk): uk:=u: vk:=v: 陰的スキームによる解 plot(x,u);

8 陽的スキームによる解 plot(x,v);

9 f1:=plot(x,u): f2:=plot(exp(-t*t),t= ,colo f3:=plot(x,uu): with (plots): display(f1,f2,f3);

10 t 非定常な問題 : 移流現象時間発展する非定常な問題として移流現象を考える 1 次元初期値問題 : 移流方程式上記の初期値問題をRiemann 問題と呼ぶがその厳密解はu0(x-at) である格子上の関数空間に対して区間 (-K,K) を分割し時間に対して区間 (0,T) を分割して各格子点上での関数値を設定する区間 (-K,K) をn 分割し区間 (0,T) をm 分割する格子点数はn+1 で

11 を定義できる格子間隔 h= 格子点の位置時刻格子点での u の値 u(,)= 陽的差分スキーム風上スキーム Lax-Wendroff スキーム Lax-Wendroff スキームの導出 Taylor 展開して u(x,t+k)=u(x,t)+k (x,t)+ ここでを代入して u(x,t+k)=u(x,t)-ak (x,t)+ これより導出される高精度スキームになっている差分スキームのプログラミングステップ 1 格子点の生成空間の区間 (-1.0,1.0) をn=200 分割 a 1 とする n:=200; gridpoints:=n; h:=evalf(2.0/gridpoints,9); deltat:=evalf(0.125*h,9); a:=1.0; (20) (21) (22) (23) (24)

12 xx:=vector(1..n+1): for i from 1 to n+1 do: xx[i]:=evalf(-1.0+h*(i-1),9); (24) ステップ 2 初期値の設定 u0:=vector(1..n+1): uexact:=vector(1..n+1): uk:=vector(1..n+1): vk:=vector(1..n+1): u:=vector(1..n+1): v:=vector(1..n+1): for i from 1 to n+1 do: if xx[i]<=0.0 then u0[i]:=1.0 else u0[i]:=0.0: end end do; if: uk:=u0: vk:=u0: m:=200; at:=evalf(a*deltat*m,9); ステップ 3 厳密解 (25) (26) for i from 1 to n+1 do: if xx[i]<=at then uexact[i]:=1.0 else uexact[i]:=0.0: end if: end do; k:=1; l:=1; ステップ4 解法 upwind (27) (28) while (k<m ) do: for i from 2 to n do: u[i]:=uk[i]-(a*deltat/h)*(uk[i]-uk u[1]:=1.0: u[n+1]:=0.0: uk:=u: k:=k+1 ステップ4 解法 Lax-Wendroff while (l<m ) do:

13 for i from 2 to n do: v[i]:=vk[i]-(0.5*a*deltat/h)*(vk[i (deltat^2)/(h^2))*(vk[i+1]-2.0*vk[i]+ v[1]:=1.0: v[n+1]:=0.0: vk:=v: l:=l+1: f1:=plot((xx,v),x= ,y= f2:=plot((xx,u),x= ,y= f3:=plot((xx,uexact),x= ,y=-0 with (plots): display(f1,f2,f3); 1 y 0 1 x

板バネの元は固定にします x[0] は常に0です : > x[0]:=t->0; (1.2) 初期値の設定をします以降 for 文処理のため空集合を生成しておきます : > init:={}: 30 番目 ( 端 ) 以外については初期高さおよび初速は全て 0 にします初期高さを x[j]

板バネの元は固定にします x[0] は常に0です : > x[0]:=t->0; (1.2) 初期値の設定をします以降 for 文処理のため空集合を生成しておきます : > init:={}: 30 番目 ( 端 ) 以外については初期高さおよび初速は全て 0 にします初期高さを x[j] 機械振動論固有振動と振動モード本事例では板バネを解析対象として数値計算 ( シミュレーション ) と固有値問題を解くことにより振動解析を行っています実際の振動は振動モードと呼ばれる特定パターンが複数組み合わされますが各振動モードによる振動に分けて解析を行うことでその現象を捉え易くすることが出来ますそこで本事例ではアニメーションを活用した解析結果の可視化も取り入れています板バネの振動