スポンサーサイト

上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。

Neumann境界条件の差分法への適用

境界での導関数の値を指定する、Neumann(ノイマン)境界条件の下での数値解析方法について考えます。

今回はDirichlet条件のときほどあっさりとはいかず、ちょっとした式変形が要ります。

Tags : 数学 数値解析

Neumann条件とは、境界条件の一種で、境界における(位置の)導関数の値を指定するタイプのものを指します。

1次元の熱伝導方程式
∂u/∂t = α・∂2u/∂x2 (α>0, 0≦x≦L) …… (*-1)
を例にとると、一般のNeumann条件は時間の関数 a(t), b(t) を用いて
∂u/∂x|x=0 = a(t) …… (*-2)
∂u/∂x|x=L = b(t) …… (*-3)
と書かれます。


(*-1)~(*-3)を差分法で解くために、適当な差分スキームを用いて式を差分方程式に変換します。

たとえば、FTCSスキームを採用した場合は
ud(xi,tn+1) = (1-2v)・ud(xi,tn) + v・{ ud(xi-1,tn) + ud(xi+1,tn) } …… (#-1)
 ( i = 1, 2, …, N-1 )

となります。ただし
v := α・⊿t/(⊿x2)

ud(xi,tn) : u(x,t)を離散化した関数 (普通、uinと書くことが多い)
⊿x : xを離散化する間隔 (サンプリング間隔)
⊿t : tを離散化する間隔 (サンプリング間隔)
xi : xを間隔⊿xで離散化したときのi番目 ( xi = i・⊿x )
tn : tを間隔⊿tで離散化したときのn番目 ( tn = n・⊿t )
N : xを分割する数( L = N・⊿x が成り立つ )

この変換は以前から何度もやっているので説明略。(参照)


さて、Neumann条件(*-2),(*-3)を差分方程式で書くとどうなるでしょうか。

単純に(*-2)式を前進差分(x-1がないから後退差分がつかえない)で近似すると
{ ud(x1,tn) - ud(0,tn) }/⊿x = a(tn)
⇔ ud(x1,tn) - ud(0,tn) = a(tn)・⊿x
⇔ ud(x0,tn) = ud(x1,tn) - a(tn)・⊿x

となります。
しかし、この方法で計算したところ、計算結果は真値から少し外れてしまいました。
つまり、この方法は、正しいかもしれないがあまり精度がよくない方法だといえます。

より正確に計算する方法を探すために、ud(x1,t) を x=x0 のまわりでTaylor展開してみます。
3次以上の項は無視し、また、見やすくするために第2引数 "tn" を省略します。
ud(x1) = ud(x0)
+ ∂ud(x)/∂x|x=x0・⊿x + (1/2)・∂2ud/∂x2|x=x0・(⊿x2)


右辺第2項は(*-2)より a(tn) に等しく、右辺第3項は(*-1)より時間微分に書き換えることができます。
ゆえに
ud(x1) = ud(x0)
+ a(tn)・⊿x + { ⊿x2/(2α) }・∂ud(x0)/∂t
⇔ ∂ud(x0)/∂t
 = { 2α/(⊿x2) }・{ ud(x1) - ud(x0) - a(tn)・⊿x }

(↑この式までは一般の差分スキームで成り立つ)
左辺を前進差分近似(FTCSだから)して
{ ud(x0,tn+1) - ud(x0,tn) }/⊿t
= { 2α/(⊿x2) }・{ ud(x1) - ud(x0) - a(tn)・⊿x }

v = α・⊿t/(⊿x2) を用いて式を整理すると
ud(x0,tn+1) = (1-2v)・ud(x0,tn) + 2v・ud(x1,tn) - 2v・a(tn)・⊿x …… (#-2)
となり、x=x0 での時間発展の式が得られました。

x=xNについても同様にして
ud(xN,tn+1) = (1-2v)・ud(xN,tn) + 2v・ud(xN-1,tn) + 2v・b(tn)・⊿x …… (#-3)
という式が得られます。

(#-1)~(#-3)までを行列でまとめると
\begin{pmatrix} 1+2v & -2v & 0 & \ldots & \ldots & 0\\ -v & 1+2v & -v & \ddots &  & \vdots\\ 0 & -v & 1+2v & -v & \ddots & \vdots\\ \vdots & \ddots & \ddots & \ddots & \ddots & 0\\ \vdots &  & \ddots & -v & 1+2v & -v\\ 0 & \ldots & \ldots & 0 & -2v & 1+2v \end{pmatrix}\begin{pmatrix} u_d\left ( x_0, t_{n+1} \right )\\ u_d\left ( x_1, t_{n+1} \right )\\ u_d\left ( x_2, t_{n+1} \right )\\ \vdotsa\\ u_d\left ( x_{N-1}, t_{n+1} \right )\\ u_d\left ( x_N, t_{n+1} \right ) \end{pmatrix}\\\ =\begin{pmatrix} u_d\left ( x_0, t_n \right )\\ u_d\left ( x_1, t_n \right )\\ u_d\left ( x_2, t_n \right )\\ \vdots\\ u_d\left ( x_{N-1}, t_n \right )\\ u_d\left ( x_N, t_n \right ) \end{pmatrix}+2vh\begin{pmatrix}-c\left ( t_n \right )\\ 0\\ 0\\ \vdots\\ 0\\ d\left ( t_n \right )\end{pmatrix}
行列の左上と右下で、行列からはみ出る分の -v が、折り返されて -2v となっています。
特に c(t)=d(t)=0 であれば、境界からの熱の流入(or流出)がなく総熱量は保存している、すなわち断熱条件を表していることがわかります。

この方法で計算してみると、先の方法より真値に近い結果が得られました。


また、完全陰解法やCrank-Nicolson法についても、同様の式変形でNeumann条件に対応した差分方程式を作ることができます。
スポンサーサイト

コメントの投稿

非公開コメント

プロフィール

fura2

Author : fura2
数学・コンピュータを中心に、考えたこと・やったことを書いていきます。

誤植等を含め、間違いはご指摘いただければ幸いです。

FC2カウンター
検索フォーム
最新記事
最新コメント
最新トラックバック
月別アーカイブ
リンク
RSSリンクの表示
上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。