スポンサーサイト

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

FTCSスキームの安定性

差分スキーム、特にFTCSスキームの安定性についての話。
実は、数値の組み合わせによってはうまく計算できないことがあるんです。

今回は行列とか固有値とかが登場します。
少なくとも固有値の求め方を知っていないと、きついかも。

ちなみに、FTCSスキームに話を限定するのは、私が一般論を知らないから
今も、わからないところは調べ調べ書いてます。

Tags : 数学 数値解析 線形代数
数式は基本的に文字で書きますが、行列だけはやむなく画像を使います。


一般に、微分方程式を差分化して数値解を計算するとき、いつも正しい解が計算できるとは限りません。
場合によっては、解が大きく振動してしまうこともありえます。
「差分スキームが安定である」とは、
計算ステップを重ねても、解が∞へ発散してしまわないことをいいます。
厳密な定義が手元の本に載ってなかったので大雑把に。
差分法で数値解析をするうえでは、用いる差分スキームの安定性に気を配る必要があります。

※ 安定だからといって、その解が正しい解とは限りません。ただ、計算結果が発散しないと言っているだけです。ややこしいですね。
※ ステップを細かくしたときに解が元の微分方程式の解に近づく性質を収束性といいます。これについては別の機会に。


差分スキームの安定性を調べる方法には大きく分けて
・行列を用いる方法
・Fourier級数を用いる方法
の2種類があります。
ここでは、私にとってわかりやすかった前者の方法を採ります。


さて、今回も
熱伝導方程式 : ∂u/∂t = α・∂2u/dx2 ( α>0 は定数 )
を例にとって進めます。
ただし、次の
初期条件 : u(x,0) = u0(x) ( u0(x) は既知の関数 )
境界条件 : u(0,t) = u(L,t) = 0 ( L>0 は定数 )
を課します。境界条件というからには、もちろん x の定義域は 0≦x≦L です。
これは、前回前々回で扱った問題を一般化した形ですね。

FTCSスキーム(時間については前進差分、空間については中央差分)による差分表示は
熱伝導方程式 :
ud(xi,tn+1) = (1-2v)・ud(xi,tn) + v・{ ud(xi-1,tn) + ud(xi+1,tn) } ……(1)
( v := α・⊿t/(⊿x2) )
初期条件 : ud(xi,0) = u0(xi) ……(2)
境界条件 : ud(0,tn) = ud(L,tn) = 0 ……(3)
となります。
udはuを離散化(≒差分化)した関数の意味です。

xi : xを離散化したときのi番目の値
tn : tを離散化したときのn番目の値
⊿x : xを離散化する間隔(サンプリング間隔)
⊿t : tを離散化する間隔

この式も以前登場しましたね。

あと、x方向のサンプリング数を M+1 としておきます。
すると、xM = M⊿x = L などが成り立ちます。


(1)式にxi = x0, x1, …, xM を具体的に代入して、行列にまとめると
\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} 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 \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}
となります。
i = 0, M のときは(3)式から x0 = xM = 0 なので、行列には入れていません。

この方程式の各部分を
t=tnでのudベクトル : ud(n)
係数行列 : A
とおくと、上式は
ud(n+1) = A ud(n)
と書けます。(行列の漸化式)
差分法で数値解析するといっても、結局は行列のかけ算をやってただけなんですね。

漸化式から一般項を求めると
ud(n+1) = A ud(n) = A2 ud(n-1) = … = An+1 ud(0)
ud(n) = An ud(0)
となり、ud(n)の値は係数行列Aのn乗と初期値ud(0)だけで決まります。


実は、ud(n) が n→∞ で発散しない(つまり、安定である)ためには、Aの固有値の絶対値がすべて1以下であればいいということがわかっています。

その理由は、Aを対角化してみるとよくわかります。
基底変換の正則行列 : P
Aの固有値を並べた対角行列 : D
とおいて、Aを対角化すると
D = P-1AP
An = PDnP-1


ここで、Aの固有値の中に絶対値が1より大のもの(λ1とする)が1つでもあれば
lim[n→∞] |λ1n| = ∞ なので
lim[n→∞] Dn (の対角要素のどれか)も発散します。
これが、解が振動してしまって不安定な状態です。

逆に、Aのすべての固有値の絶対値が1以下ならば
lim[n→∞] |λn| = 0 なので
lim[n→∞] Dn = O ( O : 零行列 )
となり、解が発散しない、安定な状態になります。
ここの議論は少し自信なし。他のサイトや文献では、差分化誤差に注目した説明が多かった...

Aが対角化できない場合の議論は知りませんが、とりあえず今回の場合は対角化できるので問題なし。

ところで、Aの固有値λは
λ = 1 - 4v・sin2( mπ/(2M) ) ( m = 1, 2, …, M-1 )
となります。
(AがM-1次なので、λもM-1個ある。)
この証明は未完成なので、完成しだい別の記事にわけて書きます。
2010/03/16追記 : 固有値の導出記事書きました。倍角公式で変形すると上の形になります。
固有値に三角関数があらわれるなんて不思議!!

このことを認めると、すべての固有値の絶対値が1以下であるためには
| 1 - 4v・sin2( mπ/(2M) ) | ≦ 1
⇔ -1 ≦ 1 - 4v・sin2( mπ/(2M) ) ≦ 1
⇔ 0 ≦ v・sin2( mπ/(2M) ) ≦ 1/2
∴ 0 ≦ v ≦ 1/2

であれば十分なことがわかります。
この条件はCFL条件(Courant-Friedrichs-Lewy条件の略)と呼ばれています。
2010/03/11追記 : CFL条件という単語は波動方程式等の特性曲線をもった一部のPDE(の差分スキーム)で定義される条件のようです。

v = α・⊿t/(⊿x2)
だったので、CFL条件安定条件は
⊿t/(⊿x2) ≦1/(2α)
つまり、⊿xを小さくしたいならその2乗のオーダーで⊿tも小さくする必要があるということです。
これは結構厳しい条件で、⊿tが小さくなると長時間たった後の状態を計算するのにとても時間がかかってしまいます。
FTCSスキームは、考え方や実装が簡単という長所がありますが、一方でこういう条件がついてしまうという短所もあるんですね。


……長かった。
次回は、これに関連して、解が発散するような条件でシミュレーションしてみます。
スポンサーサイト

コメントの投稿

非公開コメント

プロフィール

fura2

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

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

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