スポンサーサイト

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

本命 クランク-ニコルソン法 1/2

今回は、実際に広く用いられているという差分スキーム Crank-Nicolson法 を扱いました。
まとめると長くなるので記事は2分割。
この 1/2 では、Crank-Nicolson法の簡単な導出と、熱伝導方程式に適用したときの安定性について述べます。

Tags : 数学 数値解析 線形代数

以下、Crank-Nicolson法が差分スキーム(微分方程式→差分方程式の変換手法のこと)であることを強調するために、Crank-Nicolsonスキームと呼ぶことにします。

とりあえずは、次の1次元の熱伝導問題の場合に限定して考えます。
∂u/∂t = α・∂2u/∂x2 (α>0) ……(#)
初期条件 : u(x,0) = u0(x)
境界条件 : u(0,t) = u(L,0) = 0

Crank-Nicolsonスキームは、熱伝導方程式をはじめとする
∂u/∂t = ( x, t, u, ∂u/∂x, ∂2u/∂x2の関数 ) …… (*)
といった形の偏微分方程式に適用できる差分スキームです。

Crank-Nicolsonスキームは
{ ((*)式にFTCSスキームを適用したもの) + ((*)式にBTCSスキームを適用したもの) }÷2
として実現されます。

具体的に熱伝導方程式(#)にCrank-Nicolsonスキームを適用すると
{ ud(xi,tn+1) - ud(xi,tn) }/⊿t
= (α/2)・{ ud(xi+1,tn+1) - 2ud(xi,tn+1) + ud(xi-1,tn+1) }/(⊿x2)
+ (α/2)・{ ud(xi+1,tn) - 2ud(xi,tn) + ud(xi-1,tn) }/(⊿x2)
 ( i = 1, 2, …, M-1 )
となります。

ここで、各文字は
ud : uを離散化した関数 (uとは数値的に微妙に異なるので別の文字にした)
⊿x : xを離散化する間隔 (サンプリング間隔)
⊿t : tを離散化する間隔 (サンプリング間隔)
xi : xを間隔⊿xで離散化したときのi番目 ( xi = i・⊿x )
tn : tを間隔⊿tで離散化したときのn番目 ( tn = n・⊿t )
M : xを分割する数( L = M・⊿x が成り立つ )
の意味です。

v := α・⊿t/(⊿x2)として式を整理すると
(2+2v)・ud(xi,tn+1) - v・ud(xi-1,tn+1) - v・ud(xi+1,tn+1)
= (2-2v)・ud(xi,tn) + v・ud(xi-1,tn) + v・ud(xi+1,tn)
( i = 1, 2, ……, M-1 ) …… ($)

「udの各xにおける値を成分としたベクトル」(これをud(n)とおく)を用いて、($)式を行列で表示すると
\begin{pmatrix}2+2v & -v & 0 & \ldots & 0\\ v & 2+2v & -v & \ddots & \vdots\\ 0 & -v & \ddots & \ddots & 0\\ \vdots & \ddots & \ddots & \ddots & -v\\ 0 & \ldots & 0 & -v & 2+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 )\\ \vdots\\ u_d\left ( x_{M-1},t_{n+1}\right )\end{pmatrix}\\=\: \: \: \begin{pmatrix}2-2v & v & 0 & \ldots & 0\\ v & 2-2v & v & \ddots & \vdots\\ 0 & v & \ddots & \ddots & 0\\ \vdots & \ddots & \ddots & \ddots & v\\ 0 & \ldots & 0 & v & 2-2v\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 )\\ \vdots\\ u_d\left ( x_{M-1},t_n\right )\end{pmatrix}
さらに、行列
\begin{pmatrix}2-2v & v & 0 & \ldots & 0\\ v & 2-2v & v & \ddots & \vdots\\ 0 & v & \ddots & \ddots & 0\\ \vdots & \ddots & \ddots & \ddots & v\\ 0 & \ldots & 0 & v & 2-2v\end{pmatrix}^{-1}\begin{pmatrix}2+2v & -v & 0 & \ldots & 0\\ -v & 2+2v & -v & \ddots & \vdots\\ 0 & -v & \ddots & \ddots & 0\\ \vdots & \ddots & \ddots & \ddots & -v\\ 0 & \ldots & 0 & -v & 2+2v\end{pmatrix}
Aとおくと、($)式は
A ud(n+1) = ud(n) ……($$)
と書くことができます。

実際に計算するときは、この連立方程式を逐次解くことで時間ステップをすすめていきます。



[ Crank-Nicolsonスキームの安定性 ]

($$)が時間ステップを重ねても安定に計算できるか、つまりn→∞でud(n)が発散しないかを調べます。

この記事の結果と
・μ(≠0) が行列Bの固有値のとき、1/μ はB-1の固有値
・μ が行列Bの固有値, ν が行列Cの固有値のとき、μν はBCの固有値
という事実を用いると、($$)式の行列Aの固有値λは
λ = { 1 + 2v・sin2( mπ/(2M) ) }/{ 1 - 2v・sin2( mπ/(2M) ) } ( m = 1,2,…,M-1 )

ゆえに、A-1の固有値λ'は
λ' = { 1 - 2v・sin2( mπ/(2M) ) }/{ 1 + 2v・sin2( mπ/(2M) ) } ( m = 1,2,…,M-1 )

($$)式は
ud(n) = (A-1)n ud(0)
と書けるので、結局、この式がn→∞で発散しないための必要十分条件は
A-1 のすべての固有値の絶対値が1以下」となります。(参考)

そして、|λ'|≦1 はどのような v(>0), m, M についても常に成り立っているので、Crank-Nicolsonスキームは無条件で安定だということが言えます。


以上。
後半2/2では、実際のシミュレーション結果を載せ、真値との比較をやります。
スポンサーサイト

コメントの投稿

非公開コメント

プロフィール

fura2

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

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

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