スポンサーサイト

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

陰解法( 後退差分 )でもやってみる

多くの書籍の流れを汲んで、陽解法の次は陰解法でもやってみました。
今回は、完全陰解法を使った場合のシミュレーション結果とちょっとした考察とか。

Tags : 数学 物理学 数値解析 伝熱工学 線形代数
今回扱う問題の整理

 まいどおなじみ 熱伝導方程式
 ∂u/∂t = α・∂2u/∂x2 ( α>0 は定数 )
 を、差分法で解く。

 初期条件
  u(x,0) = u0(x)
 境界条件
  u(0,t) = u(1,t) = 0

 差分スキームは
 ・時間微分については後退差分スキームを
 ・空間微分については中心差分スキームを用いる。


このスキームは完全陰解法と呼ばれていますが、
前回までで扱ったFTCS(Forward Time Central Space)スキームに対応して、
BTCS(Backward Time Central Space)スキームとでも言ったほうがわかりやすいかも。
BTCSなんて言葉はほとんど使われないようだけど。

差分方程式は
{ ud(xi,tn)-ud(xi,tn-1) }/⊿t
= α・{ ud(xi+1,tn)-2ud(xi,tn)+ud(xi-1,tn) }/(⊿x2)

整理して
(1+2v)・ud(xi,tn) - v・{ ud(xi+1,tn) + ud(xi-1,tn) } = ud(xi,tn-1) ……(#)
 ( v := α・⊿t/(⊿x2) )

ud : uを差分化した関数
⊿x : xのサンプリング間隔
⊿t : tのサンプリング間隔
xi : xを離散化したときのi番目の値 ( = i⊿x )
tn : tを離散化したときのn番目の値 ( = n⊿t )

また、xをM分割したとする。( xのサンプリング数が M+1 ということ )
このとき M⊿x = 1 が成り立ちます。

(#)式を i = 1, 2, …, M-1 について書き出してまとめると
\begin{pmatrix} 1+2v & -v & 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 & \ddots & \ddots & -v\\ 0 & \ldots & \ldots & 0 & -v & 1+2v \end{pmatrix}\begin{pmatrix} u_d\left ( x_1, t_{n+1} \right )\\ u_d\left ( x_2, t_{n+1} \right )\\ u_d\left ( x_3, t_{n+1} \right )\\ u_d\left ( x_4, t_{n+1} \right )\\ \vdots\\ u_d\left ( x_{M-1}, t_{n+1} \right ) \end{pmatrix}= \begin{pmatrix} u_d\left ( x_1, t_n \right )\\ u_d\left ( x_2, t_n \right )\\ u_d\left ( x_3, t_n \right )\\ u_d\left ( x_4, t_n \right )\\ \vdots\\ u_d\left ( x_{M-1}, t_n \right ) \end{pmatrix}
Aを係数行列、ud(n)をベクトルとして
Aud(n+1) = ud(n) ……(*)

i = 0 , M のときは境界条件から常にud=0となり、この位置では明らかに安定。
ゆえに特に調べる必要もないので行列に入れていません。

BTCSスキームで計算することは、この連立方程式を解いていくことに相当します。
このような、ud(xi,tn) を計算するのに ud(xj,tn) ( j≠i ) を用いる方法を陰解法というのでした。
周りの状態をフィードバックしているんですね。

陰解法では時間ステップごとに大きなサイズの行列の計算をすることになるので、陽解法に比べて計算に時間がかかります。
その反面、陰解法のいいところは安定条件に気を使わなくてもいいということです。
結果を先に書くと、この差分スキームは無条件で安定です。

その理由は、囲みの中

 (*)式を書き直すと
  ud(n+1) = A-1ud(n)
 Aが逆行列をもつことの証明は省略。

 ところで、Aの固有値λは
  λ = 1 + 2v・sin2( mπ/(2M) )  ( m = 1, 2, …, M-1 )
 これについてはこの記事を参照。

 正則行列Aの固有値がλのとき、A-1の固有値は1/λだから
 A-1の固有値の絶対値は常に1より小

 ゆえに、BTCSスキームは無条件安定 ( 詳しい議論はこの記事を参照 )




さて、このスキームを使って計算してみた結果を載せます。

[ 諸条件 ]
 ⊿x = 0.01
 ⊿t = 0.1
 M = 100
 α = 0.01
 ud(xi,0) = 50 ( if 25≦i≦75 ), 0 ( otherwise )
 ud(0,tn) = u(1,tn) = 0

[ 結果 ]
シミュレーション結果

陽解法だと不安定になるような数値設定でも、陰解法を使えばちゃんと計算できています。
 ( α・⊿t/ )(⊿x2) = 10 > 1/2

ちなみに、連立方程式を解くアルゴリズムには LU分解 を用いました。
スポンサーサイト

コメントの投稿

非公開コメント

プロフィール

fura2

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

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

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