プロフィール

のら

Author:のら
ふわふわと日々を過ごしている学生です。
最近、人に話せる趣味が無いことに気付きました。
いい趣味募集中ですw

TOTAL

JRA 今週の重賞

競艇 重賞表示

スポーツ系ブログパーツ配布所【SBP】

最近の記事

最近のコメント

最近のトラックバック

リンク

このブログをリンクに追加する

月別アーカイブ

カテゴリー

QRコード

QRコード

ブログ内検索

RSSフィード

謙虚さが足りないんだよ。
環境変わるので、ブログ名変えました。(前:のらのなんだかなぁなブログ(仮))
スポンサーサイト
上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。

2 次元の拡散方程式の数値解析シミュレーション ( 理論編 )
さぁ、お遊びのお勉強第一弾。

偏微分方程式の数値解析!
今回は 2 次元の拡散方程式を扱います。

拡散方程式は工場の煙突のけむりのモクモクと広がり具合とか
水にインクをたらしたときの広がり具合とかを表現できるそうです。
(大学の講義面白いよ!)

拡散方程式 ( 線形 )



ただし、α は任意の定数

んで、今回は 2 次元の拡散方程式を題材にするので、
こう書き変わります(書き換えるというか同じ式だけどね)



まーこのままやってもいいんすけど、
αを 1 ってことにしちゃいますねw

気楽にやりますよwwww

---------------- 題材にする微分方程式 ------------------------------------

------------------------------------------------------------------------------

1.差分方程式の導出

これからやることを、簡単に説明すると、

・ケムリのもくもく具合を知りたい
      ↓
・でも↑の微分方程式を厳密に計算するのはめっちゃ大変
      ↓
・じゃあ点で区切って、その点のケムリの具合がわかればいいや
      ↓
・じゃあ微分もその点のケムリのぐあいの差とかで表現できんじゃね?

ってことですw

じゃあやっていきますw
まず、x y 平面を離散化して、格子状にメッシュを切ります。

そしてそれらの点を i 番目、 j 番目って名前をつけます。

んで、( i , j ) の点のケムリの具合の時間微分 (∂ u / ∂ t のことだね!)は
1個未来の同じ点のケムリのぐあいとの差を使って



ってかけます。
同様に、 x y についても差分化します。

これはちょっとコツが必要で、
1階微分を導出する際に x なら ( i - 1/2 , j ) と ( i + 1/2 , j ) の点の1階微分を求め、
それらの差を Δx でわると簡単に求めることが出来ます。

長くなるので結果だけw
(っていうか絵を描かないで説明してるからめっさわかりにくいと思うんだ・・・ごめんなさい)





これらを元の微分方程式に代入して



このような差分方程式が導出されました。
ここで、Δx = Δy として、(格子の刻みを均等にする)
式変形すると、



こんなんになります。で、Δなんちゃらのところを μ っておくと(これを拡散数っていいます)

------------------------------------ 差分方程式 -------------------------------------

--------------------------------------------------------------------------------------

1章の結末です。
↑( n+1 乗ではなくて、時刻 t での i ,j 番目の点のけむりの具合って思っていただけると良いと思います )
スタートの式と比較すると、跡形も無く見た目が変わってしまっています

以降はこの方程式の安定条件について議論します

2.安定性解析


これからの概略説明。

↑の差分方程式を計算したら値がぶっ飛んでいかないかチェックする
           ↓
ぶっ飛びそうなら式を修正する(今回はうまく行ったので、修正しないんですけどねw)

なんでこんなことをするのかと言うと、
Lax の同等定理 っていうのがありまして、

Δt とΔx とかを 0 に近づけたら元の偏微分方程式に戻る(適合性があるっていいます。)差分方程式が
安定である(←これから証明する)とき、
格子を小さく取って計算する(精度よく計算する)と元の微分方程式の解を得ることができる(収束性)

っていう定理があるんです
だからこんなしちめんどくさいことをするんですな

さぁ、うだうだ言わずにやっていきますかw
フォン・ノイマンの安定性解析って手法を用います。



(ただし i は虚数)
とおいて差分方程式に代入、整理すると
(フーリエ級数みたいにいろんな周波数の波形の和で表されたとき、どんな挙動をするか解析するって感じかな)



こんなんになります。
ここで、普通は g は複素数になるので、振幅と、位相について考えてあげないといけないのですが、
今回は実数ということで、振幅についてだけ考えます。

等比級数みたいなので習った方もいらっしゃると思いますが、
比が 1 超えると発散するみたいなのありましたよね?

アレと同じで、g の大きさが1を超えると発散すると考えます

今回は実数より、



これが条件です。
(大きさ(絶対値)で考えないといけないのでマイナスも考慮しないといけません)
この範囲でμをどこまで値を変えていいのかを考えるってことです。

ごめんね、本当にわかりにくくて

・・・まず、こっちの不等式を考えますw


この条件を満たすμは+でないといけないってことがわかり・・・ますよね?w

・・・おkだね!w



次!w


これを式変形して


すべてのθ、φについてこの条件を満たすμにしないといけないので、
右辺の最小値がμの最大値ってことになります。
こっちがめんどくさくもあり、たのしくもあり・・・

↑で言ったとおり、


の最小値を探すってことをします。
極値の条件より、(これは偏微分の教科書を参照して下さいな)







これより、極値となりえる座標値は

(0,0) , (π,0) , (0,π) , (π,π)
ってことがわかります。

ここで、ヘッセ行列式を計算するため、もっかい微分してあげます。







これを使って候補の点を判定すると、
ヘッシアンの計算は略させていただきます・・・w
知りたい人はヘッシアンで調べてみてね

(0,0) は極大
(0,π) , (π,0) は鞍点
(π,π) は極小

という結果になります。

hoge.jpg

実際に↑の関数をプロットした例がこれ。
特徴を調べたいのですが、(0,0) の極大値がでか過ぎて、他の座標値の値がどんなになってるかわかりません。

っつーことで、


って関数を定義して(↑の関数の逆数取っただけ)、
この関数 g が (π , π)で最大になってるかをチェックします。

hoge2.jpg

確かに、 (3.14 , 3.14) あたりで一番最大になってるのがわかりますねぇ
つーことで、g の最大値は (π,π)ってことがわかったんで、f の最小値も(π,π)ってことがわかりました。

ってことで、話を元に戻すと
f の最小値は 1/4 なので、拡散数 μ のとりえる値の範囲は

となります。

↑こうすれば、差分方程式が安定になるってことね

おーわり!
次の記事で、実際にプログラム組んで偏微分方程式を解いてみます
おたのしみにぃー
スポンサーサイト
テスト終了。そして・・・
やっとこさ、テストも終わりまして、夏休み突入でございます

2日連続で雀荘で徹マンを楽しみまして、
これからは切り替えて勉強に打ち込もう!ということで意気込んでおりますw

TOEIC 500 行かないとなぁ・・・
1ヶ月お酒、麻雀、競馬とかも楽しみつつ、しっかりお勉強しますねw

やっぱあれだね
自分でやろう!って思ってやる勉強(って言うか遊び?)は楽しいね

偏微分方程式いじってるのが楽しくて楽しくて(きもちわr
むしゃくしゃして買った。今も反省はしていない
H海道○幌市に住む、無職 ○○ ○○ 容疑者(20)が大量の本を買ったとして逮捕された。
調べによると、○○容疑者は、容疑を大筋で認めており、
「むしゃくしゃして買った。今も反省はしてない」と供述しているという。

ーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーー
テスト勉強の腹いせに、本買ってやった。
完全にストレス発散で買ったからなぁ・・・あとで頑張って読まないと

微分方程式による計算科学入門微分方程式による計算科学入門
(2004/02)
三井 斌友斉藤 善弘

商品詳細を見る


偏微分方程式の差分解法 (東京大学基礎工学双書)偏微分方程式の差分解法 (東京大学基礎工学双書)
(1994/06)
高見 穎郎河村 哲也

商品詳細を見る


線型代数学 (数学選書 (1))線型代数学 (数学選書 (1))
(1974/01/20)
佐武 一郎

商品詳細を見る



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