スポンサーサイト
-- / -- / -- ( -- )
上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。
--:--:-- | スポンサー広告 | page top↑
偏微分方程式が常微分方程式になる
2008 / 03 / 08 ( Sat )

二つの独立変数からなる2階の偏微分方程式

A2u/∂x2 +2B2u/∂xy+ C2u/∂y2 = f ( x, y, u, ∂u/∂x, ∂u/∂y )

は、特性方程式

A(dy/dx)2 - 2 B (dy/dx) + C = 0

と、次の方程式

A (dp/dx)(dy/dx) + Cdq/dx - f dy/dx = 0

が成り立つことと等価であることを紹介しました。

特性方程式の根は、

dy/dx = {B + B2 - AC1/2}/A

ならびに、

dy/dx = {B - B2 - AC1/2}/A

でしたから、係数A,B,C が定数なら、dy/dx も定数になります。

λ = dy/dx

とおけば、

A λ (dp/dx) + Cdq/dx - f λ = 0

となり、独立変数がx だけですので、この式は、常微分方程式になります。

結局、偏微分方程式が常微分方程式に簡単化されました。

テーマ:自然科学 - ジャンル:学問・文化・芸術

11:50:17 | 計算数理科学 | トラックバック(0) | コメント(0) | page top↑
特性方程式から型の分類
2008 / 03 / 07 ( Fri )

2つの独立変数からなる2階の偏微分方程式

A2u/∂x2 +2B2u/∂xy+ C2u/∂y2 = f ( x, y, u, ∂u/∂x, ∂u/∂y )

から、特性方程式

A(dy/dx)2 - 2 B (dy/dx) + C = 0

を導出しましたが、

これは、dy/dxの2次方程式になっています。したがって、根は、

dy/dx = {B + B2 - AC1/2}/A

ならびに、

dy/dx = {B - B2 - AC1/2}/A

になります。

ところで以前、2階偏微分方程式は、楕円型、放物型、双曲型の3種類の型に分類されることを説明しました。そのときは、

D = B2 - AC としますと、D の符号に応じて、

D > 0 なら、双曲型(Hyperbolic)

D = 0 なら、放物型(Parabolic)

D < 0 なら、楕円型(Elliptic)

となることを紹介しました。

実は、このの値は上記根の平方根の中にある値です。

すなわち、D > 0なら、二つの実根、D = 0 なら、重根、D < 0 なら二つの複素根になります。

まとめますと、

二つの独立変数からなる2階偏微分方程式は、その特性方程式が、

二つの実根を持つ場合は、双曲型

重根を持つ場合は、放物型

二つの複素根を持つ場合は、楕円型

に分類されます。

テーマ:自然科学 - ジャンル:学問・文化・芸術

21:30:16 | 計算数理科学 | トラックバック(0) | コメント(0) | page top↑
特性方程式
2008 / 03 / 06 ( Thu )

二つの独立変数からなる2階の偏微分方程式を一般形で定義すると、

A2u/∂x2 +2B2u/∂xy+ C2u/∂y2 = f ( x, y, u, ∂u/∂x, ∂u/∂y )

のように書くことができます。

いま、偏導関数を、

u/∂x = p, ∂u/∂y = q, ∂2u/∂x2 = r, ∂2u/(∂xy ) = s, ∂2u/∂y2 = t

 とおいて、上式に代入します。

すると、

Ar + 2Bs + Ct = f

になります。

ところで、p, q の全微分は、

dp = ∂p/∂x ・ dx + p/∂y ・ dy = r dx + s dy

dq = ∂q/∂x ・ dx + q/∂y ・ dy = s dx + t dy

で定義されますので、これらから、r t を逆算して上式に代入すると、

A(dp - sdy) /dx + 2 B s + C(dq - sdx) = f

が得られます。これをさらに s で整理すると、

s{A(dy/dx)2 - 2 B (dy/dx) + C} - { A (dp/dx)(dy/dx) + Cdq/dx - f dy/dx } = 0

となります。

この式がすべての s に対して常に成り立つためには、

A(dy/dx)2 - 2 B (dy/dx) + C = 0

でかつ、

A (dp/dx)(dy/dx) + Cdq/dx - f dy/dx = 0

でなければなりません。

前者の式のことを、特性方程式(Characteristics Equation)と呼びます。

テーマ:自然科学 - ジャンル:学問・文化・芸術

21:37:55 | 計算数理科学 | トラックバック(0) | コメント(0) | page top↑
三次元熱伝導方程式をSOR法で解く(その3)
2008 / 03 / 05 ( Wed )
計算例として、立方体内外の熱伝導を計算するFortranプログラムを紹介します。
立方体内外の熱伝導率はそれぞれ10、
1 と設定します。
境界条件は、入口境界i =1 面 でT =100度、出口境界i = IF 面でT =50度 とします。
それ以外の外部境界は、法線方向に
T/∂n=0 となるようノイマン境界条件を与えす。
初期値として、入口境界以外の点で
T =50度を与えます。
格子間隔は
x =y = z = 0.1で、過緩和係数はω = 1.5とします。

これら境界条件・初期条件に設定したFortranプログラムは、下記のページに置いてありますので、ご自由のお使いください。


三次元熱伝導方程式のFotranプログラム

テーマ:自然科学 - ジャンル:学問・文化・芸術

21:31:46 | 計算数理科学 | トラックバック(0) | コメント(0) | page top↑
三次元熱伝導方程式をSOR法で解く(その2)
2008 / 03 / 04 ( Tue )
クランク・ニコルソン法で差分近似された式に、SOR法を適用すると、

L( Tn+1i,j,k)m+1=(1 - ω ) L(Tni,j,k)m + ω κ [x {(Tn+1i+1,j,k)m(Tn+1i-1,j,k) m+1}/2 + y{(Tn+1i,j+1,k)m (Tn+1i,j-1,k )m+1}/2+ z{(Tn+1i,j,k+1)mTn+1i,j,k-1 )m+1}/2+ di,j,k]

式が複雑になりわかりにくいですが、

di,j,k = Tni,j,k + x (Tni+1,j,k - 2Tni,j,kTni-1,j,k) 2 + y(Tni,j+1,k - 2Tni,j,k Tni,j-1,k )/2+ z(Tni,j,k+1 - 2Tni,j,k Tni,j,k-1 )/2

また、

x = ∆t/(∆x)2, ℓy = ∆t/(∆y)2, ℓz = ∆t/(∆z)2

L=1+x+y+z

式をさらに変形すると、

( Tn+1i,j,k)m+1=(Tni,j,k)m + ω κ [x {(Tn+1i+1,j,k)m(Tn+1i-1,j,k) m+1}/2L + y{(Tn+1i,j+1,k)m (Tn+1i,j-1,k)m+1}/2L+ z{(Tn+1i,j,k+1)m(Tn+1i,j,k-1)m+1}/2L+ di,j,k/L - (Tni,j,k)m]

となり、SOR法で反復計算しながらTn+1i,j,k を求める陰解法の差分計算式の形ができます。

テーマ:自然科学 - ジャンル:学問・文化・芸術

21:37:44 | 計算数理科学 | トラックバック(0) | コメント(0) | page top↑
三次元熱伝導方程式をSOR法で解く(その1)
2008 / 03 / 03 ( Mon )
2次元熱伝導方程式をExcelで解く方法を紹介しましたが、3次元になると、さすがにExcelで計算するのは現実的ではありません。

ここでは、三次元熱伝導方程式をクランクに・ニコルソン法で差分近似して、SOR法で解く、Fortranプログラムを紹介します。

まず、三次元熱伝導方程式を次式で定義します。

T/∂t = κ(2T/∂x2 + ∂2T/∂y2+ 2T/∂y2)

T は温度です。これをクランク・ニコルソン法で差分近似すると、

( Tn+1i,j,kTni,j,k)/t = κ {(Tn+1i+1,j,k - 2Tn+1i,j,kTn+1i-1,j,k) / (∆x)2 + (Tn+1i,j+1,k - 2Tn+1i,j,k Tn+1i,j-1,k ) / (∆y)2 + (Tn+1i,j,k+1 - 2Tn+1i,j,k Tn+1i,j,k-1 ) / (∆z)2 +(Tni+1,j,k - 2Tni,j,kTni-1,j,k) / (∆x)2 + (Tni,j+1,k - 2Tni,j,k Tni,j-1,k ) / (∆y)2 + (Tni,j,k+1 - 2Tni,j,k Tn+i,j,k-1 ) / (∆z)2 }/2

とたいへん長い式になります。

テーマ:自然科学 - ジャンル:学問・文化・芸術

21:35:03 | 計算数理科学 | トラックバック(0) | コメント(0) | page top↑
Excelで2次元熱伝導方程式を解く(その4)
2008 / 03 / 02 ( Sun )

マクロプログラムを実行した結果を示します。

このマクロプログラムでは、40時間ステップ計算するように設定されており、その間、領域3の2次元セル領域C21:M31の値を、領域2の2次元セル領域C7:M17へのコピーを繰り返し、かつセル番号E4にあるTimeにセル番号B4にあるDelta tの値を加算します。

Timeの値は最終的には、Time=0.25x40=10となります。時間ステップの数を増やしたい場合には、繰り返し回数の数値を変更するか、繰り返しマクロプログラムを実行します。  

40時間ステップ計算したときの計算結果は次のようになりました。

Excelsheet2.png

領域1のTimeの値が10に更新され、領域2,3のセルにも計算結果が表示されています。さらに、計算された結果は、領域4に可視化されます。

さらに、Timeが40になるまで繰り返し計算した場合の結果をも示します。80時間ステップでほぼ収束解(定常解)が得られました。

Excelsheet3.png

ここで紹介したExcelのプログラムは、下記ページに置いてありますので、ご自由にお使いください。

二次元熱伝導方程式のExcelプログラム

テーマ:自然科学 - ジャンル:学問・文化・芸術

07:52:59 | 計算数理科学 | トラックバック(0) | コメント(0) | page top↑
Excelで2次元熱伝導方程式を解く(その3)
2008 / 03 / 01 ( Sat )

ここで紹介する方法で使用する簡単なマクロプログラムを以下に示します。

Sub 二次元熱伝導方程式()
For i = 1 To 40
Range("C21:M31").Copy
Range("C7:M17").PasteSpecial (xlPasteValues)
Range("E4").Value = Range("E4").Value + Range("B4").Value
Next i
End Sub

このプログラムでは、
領域3から領域2へのコピーと、時間ステップの更新が行われます。

このプログラムを、Excelのツールバーにある、ツール/マクロ/Visual Basic Editorを立ち上げて、左ボックスに表示された、Sheet1をダブルクリックして開き、その中にコピーします。

マクロプログラム

最終的に、ツール/マクロ/マクロで、作成済みのマクロプログラム名を選択して、実行ボタンを押せばマクロプログラムが実行されて計算が開始します。

テーマ:自然科学 - ジャンル:学問・文化・芸術

11:09:20 | 計算数理科学 | トラックバック(0) | コメント(0) | page top↑
Excelで2次元熱伝導方程式を解く(その2)
2008 / 02 / 29 ( Fri )

Excelのシートを最大限活用します。ここで紹介する手法は、すでにインターネット上では、いくつか紹介されているものです。


まず、今回用意したExcelシートを、次に示します。

Excelsheet1.png

シートは大きく分けて、初期値を与える「領域1」、n 時間ステップの計算結果を保管する「領域2」、n+1時間ステップの計算ならびに計算結果を保管する「領域3」、そして計算結果を可視化する「領域4」の4つの領域からなります。

領域1では4つの変数、∆t = 0.25、x = 0.1、κ = 0.01そして初期時間、 t = 0を規定しています。

次に、領域2には、計算格子の各点におけるn 時間ステップの値が表示されます。ただし、ここでは計算格子点数を、IF = 11、JF = 11としました。計算初期値として、セル番号C7からC17まで1を代入し、その他のセルにはすべて0が代入されています。

領 域3がn+1時間ステップの計算をしている中枢部分であり、境界以外の内点のセルで、差分近似式が解かれます。

たとえば、セル番号D22をクリックすると、

fx=D8+$B$4*$D$4*(E8+C8+D9+D7-D8*4)/$C$4/$C$4

という演算式が現れます。Excelでセルの行番号と列番号にそれぞれ$マークが付いた場合には、固定されたセル番号を参照しますので、$B$4、$C$4、$D $4はそれぞれ領域1のDelta t、Delta x、Kappaの値である0.25、0.1、0.01を参照します。
また、セル番号D8、E8、C8、D9、D7は、それぞれ、
uni,j , uni+1.j , uni-1,j , uni,j+1  ,  uni,j-1 の値が保 管されている領域2のセル番号に相当します。
したがって、たとえばセル番号D22においては、格子点 i = 2、j = 2における値が計算されます。

一方、境界では、ノイマン境界条件が与えられている上辺と底辺に相当する格子点において計算が必要になります。

ここでは、計 算する差分近似された式は陽解法ですので、n+1 時間ステップで計算された1つ内点の値を用いて外挿します。たとえば、上辺 境界に位置するセル番号D21の計算式は、fx= D22と置きます。

境界と内点の計算式を用いた計算がすべて終了した段階で、領域3にはn+1時間ステップの計算結果 が保管されます。

テーマ:自然科学 - ジャンル:学問・文化・芸術

21:54:32 | 計算数理科学 | トラックバック(0) | コメント(0) | page top↑
Excelで2次元熱伝導方程式を解く(その1)
2008 / 02 / 28 ( Thu )
今回から何回かに分けて、2次元熱伝導方程式を Excel を用いて解く方法について紹介します。

1次元熱伝導方程式については、すでにExcelを用いて計算した例を紹介しました。1次元の場合には、画面上でExcelを手動で操作しながら簡単に計算することができます。しかしながら、2次元の場合には無理かもしれません。

そこで、簡単なマクロプログラムを追加するだけで計算できる方法について紹介します。ここでの内容は、すでに、ある学会から依頼されて執筆し掲載された原稿に基づいています。

まず、2次元熱伝導方程式を次のように定義します。

u/∂t = κ(2u/∂x2 + ∂2u/∂y2)

これを陽解法で差分近似すると、

( un+1i,juni,j)/t = κ {(ui+1.j - 2ui,jui-1,j) / (∆x)2 + (ui,j+1 - 2ui,j ui,j-1 ) / (∆y)2 }

簡単にするため、∆x = ∆y =κt/ (x )2 とおけば、

un+1i,j = uni,j + (uni+1.juni-1,j + uni,j+1uni,j-4uni,j)

のように比較的簡単な式を導出することができます。

具体的な熱伝導問題として、次の図に示す境界条件を設定します。

境界条件

クランク・ニコルソン法を、ガウス・ザイデル法やSOR法により数値計算する手法を紹介しましたが、マクロプログラムにすると結構複雑になり、Excel 本来の利点が生かせないので、ここでは最も単純な陽解法による方法を紹介していきます。

テーマ:自然科学 - ジャンル:学問・文化・芸術

21:31:08 | 計算数理科学 | トラックバック(0) | コメント(0) | page top↑
前ページ | ホーム | 次ページ
上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。