スポンサーサイト
-- / -- / -- ( -- )
上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。
--:--:-- | スポンサー広告 | page top↑
双曲型方程式を特性曲線法で具体的に解く(その3)
2008 / 03 / 17 ( Mon )

次に、点NにおけるpN, qNの値を求めます。

その際に、

AL λL (pN-pL)/(xN-xL) + C L(qN-qL)/(xN-xL)= 0

AM λM (pN-pM)/(xN-xM) + CM(qN-qM)/(xN-xM)= 0

を連立して解く必要があります。ただし、これらの式は、

AL λL (pN-pL)+ C L(qN-qL)= 0

AM λM (pN-pM)+ CM(qN-qM)= 0

のように簡略化できます。

AL=AM=1

CL= - uL2 = -0.42 = -0.16, CM= - uM2 = -(-0.65)2=0.4225

pL=10*0.2=2, pM=10*0.3=3, qL=3*0.2=0.6, qM=3*0.3=0.9

ですので、結局、

0.4(pN - 2.0)- 0.16(qN - 0.6)= 0

-0.65(pN - 3.0) - 0.4225(qN - 0.9)= 0

から、pN, qNを求めれば、

pN=2.45524,  qN=1.73810

となります。

求められた、xN, yN, pN, qNを、

uN - uL = (pN+pL)(xN - xL) + (qN+qL)(yN - yL)

に代入すれば、

uN = 0.4 +  (2.0 + 2.45524)(0.26190 - 0.2) /2+ (0.6 + 1.73810)(0.024762- 0)

で、これを計算すれば、

 uN = 0.56684

となります。

これが、特性曲線法により求められた解です。

スポンサーサイト

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

10:37:22 | 計算数理科学 | トラックバック(0) | コメント(0) | page top↑
双曲型方程式を特性曲線法で具体的に解く(その2)
2008 / 03 / 16 ( Sun )

まず、偏導関数p,qは、

p =du/dx = 10x ,  q =du/dy = 3x

と求まります。次に特性方程式は、

dy/dx 2 - u2 = 0

になります。したがって、この根は、

λ =  dy/dx = u , ならびに -u

です。

つぎに、xN, yNを求めます。

λL = uL= 0.2 + 5*0.22 = 0.4

λM= - uM=-0.2 - 5*0.32 = -0.65

ですので、

yN- 0 =0.4(xN- 0.2)

yN-  0 = - 0.65(xN- 0.3)

となり、これらを連立して解けば、

xN=0.26190,  yN=0.024562

と点Nの座標が求まります。

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

13:48:17 | 計算数理科学 | トラックバック(0) | コメント(0) | page top↑
双曲型方程式を特性曲線法で具体的に解く(その1)
2008 / 03 / 15 ( Sat )

双曲型方程式

2u/∂x2 - u22u/∂y2= 0

を特性曲線法で具体的に解いてみます。

まず、

D = B2 - AC = u2 > 0

ですので、この方程式は双曲型です。

いま、x=0.2にある点Lx=0.3にある点Mの間にあり、y > 0にある特性曲線が交叉した点Nの座標 (xN,yN) におけるuN を求めてみます。

ただし、u は初期値 y=0 上で、0 x 1 のとき、

u = 0.2 + 5x2

で、かつ

u/∂y = 3x

を満足するとします。

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

16:50:40 | 計算数理科学 | トラックバック(0) | コメント(0) | page top↑
2階偏微分方程式を特性曲線法で解く(その4)
2008 / 03 / 14 ( Fri )
(pN,qN)の値が既知となりました。
もともと、p, qは、

p = u/∂x , q =u/∂

でした。
最後に、点Nにおける
uNの値を求めるために、uの全微分式、

du =
u/∂x・dx+u/∂y・dy

を差分近似します。すなわち、

uN - uL = {(pN+pL)(xN - xL) + (qN+qL)(yN - yL) }/2

ただし、
p, qの値は、点Lと点Nにおける値を平均します。
また、直線MNで成り立つ差分近似式

u
N -
uM = {(pN+pM)(xN - xM) + (qN+qM)(yN -  yM)}/2

を解いてもかまいません。解は同じになります。
これより、
点Nにおける解uNを求めることができます。
特性曲線法では、点Nの解を求める計算を他の任意の点でも行い、求められた複数の点における解を用いて、さらに次の計算を同様に繰り返すことにより、解の伝播を手計算で求めることができます。



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

23:00:18 | 計算数理科学 | トラックバック(0) | コメント(0) | page top↑
2階偏微分方程式を特性曲線法で解く(その3)
2008 / 03 / 13 ( Thu )
点Nの座標 (xN,yN) が既知となりましたが、
次に、
λLでは、

AL λL (dp/dx)L + C L(dq/dx)L - fL λL = 0

λMでは、

AM λM (dp/dx)M + C M(dq/dx)M - fMλM = 0

の各式が成り立ちますので、それぞれを、直線LN、LM上で差分近似します。

すると、

AL λL (pN-pL)/(xN-xL) + C L(qN-qL)/(xN-xL) - fL λL = 0

AM λM(pN-pM)/(xN-xM) + C L(qN-qM)/(xN-xM) - fM λM = 0

の2式が得られます。

これらの式における未知変数は、(pN,qN)のみです。

したがって、上記2式を連立することで(pN,qN)を求めることができます。

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

21:46:42 | 計算数理科学 | トラックバック(0) | コメント(0) | page top↑
2階偏微分方程式を特性曲線法で解く(その2)
2008 / 03 / 12 ( Wed )
λLでは、

λL = (dy/dx)L

ですから、これが点Nを通るとすれば、次の差分近似式が成り立ちます。

λL = (yN- yL)/(xN- xL)

少し書き換えれば、

yN- yL=λL(xN- xL)

同様に、 λMでは、

λM=  (dy/dx)M

ですから、これが点Nを通るとすれば、

λM = (yN- yM)/(xN- xM)

で、

yN- yM=λM(xN- xM)

です。

ところで、点L,Nにおける(xL,yL), (xM,yM)は既知の値です。

未知の値は、(xN,yN)です。

したがって、

yN- yL=λL(xN- xL)

yN- yM=λM(xN- xM)

の2式を連立して計算することにより、(xN,yN)を求めることができます。

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

20:48:18 | 計算数理科学 | トラックバック(0) | コメント(0) | page top↑
2階偏微分方程式を特性曲線法で解く(その1)
2008 / 03 / 11 ( Tue )
今日から何回かに分けて、特性曲線法(Method of Characteristic Curves)について紹介します。
特性曲線法とは、偏微分方程式を特性曲線上で解く方法です。
コンピュータが進歩した今日では、あまり使われない方法ですが、むしろコンピュータがまだなかった時代に、手計算やタイガー計算機を用いて、広く用いられていた方法です。特に双曲型偏微分方程式の特徴を理解する上で、大変役に立つ方法です。
準備として、次の図を用意します。

特性曲線法

図中には、2本の特性曲線が示されています。
まず、点Lから延びた
dy/dx = {B + B2 - AC1/2}/A を勾配に持った特性曲線をλL、点Mから延びた dy/dx = {B - B2 - AC1/2}/A を勾配に持った特性曲線をλMとします。さらに、これらの特性曲線が交差した点を、点Nとします。それぞれの特性曲線上では、図中に示したそれぞれ2つの式が成り立ちます。すなわち、
λLでは、

λL = (dy/dx)L

ならびに

AL λL (dp/dx)L + C L(dq/dx)L - fL λL = 0

λMでは、

λM = (dy/dx)M

ならびに

AM λM (dp/dx)M + C M(dq/dx)M - fMλM = 0

が成り立ちます。

これら4つの式を連立して解くことにより、解を求めることができます。

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

21:21:37 | 計算数理科学 | トラックバック(0) | コメント(0) | page top↑
特性曲線上の常微分方程式
2008 / 03 / 10 ( Mon )
特性曲線の図を簡略化して、改めて示します。

特性曲線その2

初期値上の2点を、それぞれL、Mとして、それぞれから伸びている特性方程式の根を勾配に持つ直線を図示します。
たとえば、L点から伝播した解を示す2つの直線の内、dy/dx = {B + B2 - AC1/2}/Aの勾配を、λL+ とします。これを言い換えれば、「特性曲線の勾配」です。
このような勾配を持つ特性曲線は、初期値のあらゆる点から伸びており、全空間を埋め尽くしています。
特性方程式の根とは、すなわち、特性曲線の勾配です。
xy
空間は、特性方程式の根を勾配に持つ特性曲線で埋め尽くされています。

ここからが重要ですが、そのような勾配をもった特性曲線上では、先に示した常微分方程式が成り立ちます。
たとえば、
λL+を勾配に持つ特性曲線上では、

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

が成り立ちます。
ここで紹介している偏微分方程式は、
その特性方程式の根を勾配に持つ特性曲線上で成り立つ常微分方程式を解くことで解を求めることができます。
偏微分方程式は、事実上、数値計算でしか解くことができませんが、常微分方程式でしたらば、解析的に解くこともできます。


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

21:23:24 | 計算数理科学 | トラックバック(0) | コメント(0) | page top↑
特性方程式の根を図示する
2008 / 03 / 09 ( Sun )
特性方程式の根をわかりやすく説明するために、下にように図示してみます。

特性曲線

独立変数 (x,y) からなる2階偏微分方程式は、xy空間上で図示することができます。
初期条件として、y = 0(y > 0) とすれば、
y = 0 上の1本の直線が初期値になります。
この直線上の任意の点から、
y > 0方向に解は伝播します。
いま、
(x1,0) (x2,0) の2点を考えます。そこから、dy/dxが一定値になる方向に直線をひくことができます。
仮にその直線の
dy/dx が特性方程式の根であるとすれば、図中のようにそれぞれの点から2本ずつ直線をひくことができます。
このような直線のことを、特性曲線(Characteristic Curve)と呼びます。
直線なのに曲線というのは変ですが、ここでの説明はあくまで、1次精度を仮定していますので、2次精度になれば直線は曲線になります。

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

10:15:33 | 計算数理科学 | トラックバック(0) | コメント(0) | 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↑
| ホーム |
上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。