1 / 52

８．４ 偏微分方程式の数値解法

８．４ 偏微分方程式の数値解法. のとき. ［偏微分の復習（１）］. 座標系の座標値とすると，. この. と. を. と. お互いに 直交する値 であるから，. は 独立 。. は 意味がない 。. そこで. お互いに 相手の変数を定数 とみなして微分する。. 偏微分 と呼ぶ。. ［記号］. の替わりに. を用いる。. 偏微分の例. を定数とみなして. を定数とみなして. を見慣れないからといって恐れないこと。 常微分より簡単 だよ。. ［常微分と偏微分の関係］. ［偏微分の復習（２）］. ［媒介変数表示］. に関して偏微分できるとき，.

aleda
Download Presentation

８．４ 偏微分方程式の数値解法

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. ８．４ 偏微分方程式の数値解法 のとき ［偏微分の復習（１）］ 座標系の座標値とすると， この と を と お互いに直交する値であるから， は独立。 は意味がない。

  2. そこで お互いに相手の変数を定数とみなして微分する。 偏微分と呼ぶ。 ［記号］ の替わりに を用いる。

  3. 偏微分の例 を定数とみなして を定数とみなして を見慣れないからといって恐れないこと。常微分より簡単だよ。

  4. ［常微分と偏微分の関係］ ［偏微分の復習（２）］ ［媒介変数表示］ に関して偏微分できるとき， が，ともに 以上の関係は，いたるところで出てくるので重要！！

  5. ［勾配］ ［ベクトル解析の復習①］ は，それぞれ 方向の単位ベクトル という演算子（ナブラ）を導入すると

  6. ２次元の地形図を思い浮かべるとよい。 ２次元の▽（ナブラ）は スカラー値 を標高値とすると， は傾斜方向ベクトルである。 イメージをつかむために 絶対値： 方向：

  7. ［発散］ベクトル場 に対して 簡単のために２次元の数値地図を考え， ［ベクトル解析の復習③］ 左のメッシュでは 不変とする と， が正のとき，傾斜方向は広が ベクトル場 がスカラー の勾配のとき り，負のときは狭まる形になる。 ここで， とあらわし これを のラプラシアンという。 （正のとき凹型地形，負のとき凸型地形） （あるいは，隣り合うベクトルの広がり具合）

  8. ［回転］ベクトル場 に対して ［ベクトル解析の復習④］ なお， と表記されることもある。

  9. 公式を丸暗記しないこと。意味を理解する。できなければ，どの資料に書いてあるかを記憶し，使うときに参照する。使っているうちに慣れてくる。公式を丸暗記しないこと。意味を理解する。できなければ，どの資料に書いてあるかを記憶し，使うときに参照する。使っているうちに慣れてくる。 ［公式］ ［ベクトル解析の復習⑤］

  10. ［ベクトル解析の復習⑥］

  11. ［ベクトル解析の復習⑦］

  12. ８．４．１　数値解法における差分の方法（１）差分の種類８．４．１　数値解法における差分の方法（１）差分の種類 常微分の場合と同様，以下のような方法がある。 　　①前進差分（forward difference） 　　②中心差分（centered difference） 　　③後退差分（backward difference）

  13. （２）前進差分 微小格子を考え，以下のように近似する。

  14. （３）後退差分 微小格子を考え，以下のように近似する。

  15. （４）中心差分 微小格子を考え，以下のように近似する。 中心差分は，前進差分から後退差分を差し引いて，１／２にした形である。

  16. （５）差分式使用上の注意点 どの方法を選択するかは， 　　①与えられた条件， 　　②問題の性格， 　　③モデル化のしやすさ 等に応じて決めること。 流体力学の分野では，空間成分に関して，速度成分が正であれば後退差分を，負であれば前進差分を使用する風上差分も用いられている。

  17. ８．４．２　物理現象と偏微分方程式（１）考え方８．４．２　物理現象と偏微分方程式（１）考え方 ①偏微分方程式は，物理現象を支配する法則等により体積無限小の極限で導出される。 ②有限な体積要素における法則により，基礎式が直接導出されることも多い。 　（例：コントロールボリューム法，直接差分法）

  18. （２）２階偏微分方程式の分類 ［分類］ ２階偏微分方程式の一般式 ① のとき楕円型 ［例］ （ラプラス方程式：Laplace’s equation） （ポアソン方程式：Poisson’s equation） ② のとき放物型 ［例］ （拡散方程式など） ③ のとき双曲型 ［例］ （波動伝播など）

  19. （ラプラス方程式） （３）楕円型のプログラム例 （ポアソン方程式） ［シートの定義］

  20. VBAによる定義　①データ宣言 Private ポアソン As Boolean Private KX As Integer Private KY As Integer Private NN As Integer Private EPS As Double Private U() As Double Private ID1 As Integer Private ID2 As Integer Private F() As Double Private UMax As Double Private UMin As Double Private UDX As Double

  21. ②初期条件の設定と境界条件設定 Private Sub 初期条件(KX, KY, MX, MY, DX, DY) For j = 0 To KY For i = 0 To KX U(0, i, j) = 0: U(1, i, j) = 0: F(i, j) = 0 If ポアソン Then ‘　ポアソンの方程式のとき右辺設定 F(i, j) = 80 * (DX * (MX - i) + DY * (MY - j)) End if Next Next End Sub Private Sub 境界条件(KX, KY, DX, DY) For k = 0 To 1 For i = 0 To KX: U(k, i, 0) = 0: U(k, i, KY) = 16: Next For j = 0 To KY: U(k, 0, j) = 0: U(k, KX, j) = 8: Next Next End Sub

  22. ③計算本体（その１） Private Sub 計算(KX, KY, NN, EPS) Dim MX As Integer: MX = KX + 1 Dim MY As Integer: MY = KY + 1 Dim DX As Double: DX = 1# / KX Dim DY As Double: DY = 1# / KY ReDim U(1, KX, KY), F(KX, KY) Dim F1 As Double: F1 = 1 / (DX * DX) Dim F2 As Double: F2 = 1 / (DY * DY) Dim F3 As Double: F3 = 0.5 / (F1 + F2) 初期条件 KX, KY, MX, MY, DX, DY 境界条件 KX, KY, DX, DY

  23. ④計算本体（その２） ' 収斂計算 ID1 = 1: ID2 = 0 For N = 0 To NN - 1 ID = ID1: ID1 = ID2: ID2 = ID: ER = 0 For j = 1 To KY - 1 For i = 1 To KX - 1 U(ID2, i, j) = F3 * (F1 * (U(ID1, i + 1, j) + U(ID1, i, j + 1)) + _ F2 * (U(ID1, i - 1, j) + U(ID1, i, j - 1)) + F(i, j)) If Abs(U(ID2, i, j)) > EPS Then E = Abs((U(ID2, i, j) - U(ID1, i, j)) / U(ID2, i, j)) If E > ER Then ER = E End If Next Next

  24. ⑤計算本体（その３） If ER < EPS Then Exit For If (N Mod 50) = 0 Then ‘計算途中経過表示 With Worksheets("データ") .Cells(2, 5) = N .Cells(2, 6) = Format(ER, "#0.000000") End With Application.ScreenUpdating = True Application.ScreenUpdating = False End If Next Application.ScreenUpdating = True End Sub

  25. ⑥データの設定 Sub データ設定() With Worksheets("データ") KX = Val(.Cells(2, 1)) KY = Val(.Cells(2, 2)) NN = Val(.Cells(2, 3)) EPS = Val(.Cells(2, 4)) End With ポアソン = False End Sub

  26. ⑦結果の設定 Sub 結果設定() With Worksheets("結果") UMax = U(ID2, 1, 1): UMin = U(ID2, 1, 1) For j = 1 To KY - 1 For i = 1 To KX - 1 .Cells(j, i) = U(ID2, i, j) If UMax < U(ID2, i, j) Then UMax = U(ID2, i, j) If UMin > U(ID2, i, j) Then UMin = U(ID2, i, j) Next Next UDX = UMax - UMin End With End Sub

  27. ⑧ボタンのClickイベントハンドラ Sub ボタン1_Click() ‘ラブラスの方程式のためのClickイベントハンドラ データ設定 ポアソン = False 計算 KX, KY, NN, EPS 結果設定 End Sub Sub ボタン3_Click() ‘ポアソンの方程式のためのClickイベントハンドラ データ設定 ポアソン = True 計算 KX, KY, NN, EPS 結果設定 End Sub

  28. ⑨結果を分布図としてセルの色で表示する処理⑨結果を分布図としてセルの色で表示する処理 Sub ボタン2_Click() Dim CD(10) As Integer CD(0) = 5: CD(1) = 41: CD(2) = 33: CD(3) = 34: CD(4) = 36 CD(5) = 6: CD(6) = 44: CD(7) = 45: CD(8) = 46: CD(9) = 3 Worksheets("分布").Select MY = KY + 1 With Worksheets("結果") For j = 1 To KY - 1 For i = 1 To KX - 1 ID = 9 * (Val(.Cells(KY - j, i) - UMin) / UDX) If ID > 9 Then ID = 9 Worksheets("分布").Cells(j, i).Select Selection.Interior.ColorIndex = CD(ID) Selection.Interior.Pattern = xlSolid Next Next End With End Sub

  29. ⑩結果分布図 ポアソンの方程式の結果 (Umax=14.93, Umin=0.03) ラプラスの方程式の結果 (Umax=15.35, Umin=0.00) ラプラスの方程式の結果の方が分布として 滑らかであることに注意されたい。

  30. （４）電磁場の偏微分方程式 （ガウスの法則）　　　　　　　　　　　ポアソンの方程式 電磁場における支配方程式 （マクスウェルの方程式） （電磁誘導の法則）　　　　　　　　 （磁束密度湧き出しなしの法則）　ラプラスの方程式 （アンペールの法則） ： 電束密度 基本式が同じでも， 条件によって 偏微分方程式が異なることを ここでは 確認する。 ： 電場の強さ ： 磁場の強さ ： 磁束密度 ： 電流密度 ： 電荷

  31. （４）電磁場の偏微分方程式 ，透磁率 ，電荷 を用いた電束密度と電場強さ，磁束密度と 磁場の強さ，電流密度と電場の強さの関係 誘電率 これらをマックスウェルの方程式に代入して 右の式を時間で偏微分して，左の式を代入すると

  32. 条件で変わる偏微分方程式 として整理すると 真空中では 同じ方程式でも， 条件によって偏微分方程式が異なる （双曲型：真空中の電場伝播） 電場の強さの 時間変化が小さいとき 定常状態のとき元の式 電位 を用いて （放物線型：拡散方程式） （楕円型：ポアソン方程式）

  33. （５）波動方程式 （梁の振動：１次元） １次元（平面波） ２次元（円柱波） ３次元（球面波） 極座標による２次元波動方程式（太鼓など円形の膜の振動のとき）

  34. １次元波動方程式の解釈 １次元（平面波） 隣り合う両節点の変位差を変位とみなせば， 加速度が変位に比例するバネと同じ

  35. プログラム例 とおいて，以下の差分方程式としてプログラミングする。 ［シートの定義］（以下のように３種類のシートとコマンドボタンを用意する）

  36. VBAによる定義　①データ宣言と初期値設定 Private Z(2, 100) As Double Private V(2, 100) As Double Private ID1 As Integer Private ID2 As Integer Sub 初期値() For i = 0 To 50 A = Sin(3.14159265359 * i / 100) Z(0, i) = A： Z(0, 100 - i) = A ‘　端点での誤差を少なくするために V(0, 100 - i) = 0: V(0, i) = 0 ‘このように設定している Next With Worksheets("結果") For i = 0 To 100 .Cells(1, i + 2) = i .Cells(2, i + 2) = Z(0, i) Next End With End Sub

  37. VBAによる定義　②処理本体（その１） Sub 波動方程式1次元(NumLoop, NN, SaveN, DT, DX, Alfa) Beta = Alfa * Alfa / (DX * DX): ID1 = 0: ID2 = 1: KK = 0 For k = 1 To NumLoop For j = 1 To NN - 1 Z(ID2, j) = Z(ID1, j) + DT * V(ID1, j) Acc = Beta * (Z(ID1, j + 1) - 2 * Z(ID1, j) + Z(ID1, j - 1)) V(ID2, j) = V(ID1, j) + DT * Acc Next ‘　境界条件（両端を固定とする） V(ID2, 0) = 0: V(ID2, 100) = 0 Z(ID2, 0) = 0: Z(ID2, 100) = 0

  38. VBAによる定義　③処理本体（その２）とボタンのClickイベントハンドラVBAによる定義　③処理本体（その２）とボタンのClickイベントハンドラ If (k Mod SaveN) = 0 Then ‘ 指定された計算回数間隔で途中結果を保存する KK = KK + 1 With Worksheets("結果") .Cells(KK + 2, 1) = KK * SaveN * DT For j = 0 To NN .Cells(KK + 2, j + 2) = Z(ID2, j) Next End With End If Temp = ID2: ID1 = ID2: ID2 = Temp Next End Sub Sub ボタン2_Click() 初期値 波動方程式1次元 40000, 100, 100, 0.01, 0.05, 0.05 End Sub

  39. ⑨結果グラフ 両端での速度・位置を固定しているので， 振幅が次第に小さくなる。 中央の振幅 T=100までの弦の変化 αの値を変化させて， 色々確かめてみよう。 また，αが大きすぎると 振動状態になってしまうことも 確認すること。

  40. （６）保存量の偏微分方程式 境界面を移動する流束ベクトル 熱拡散方程式 の条件による定常状態では

  41. プログラム例 ［シートの定義］（以下のように３種類のシートとコマンドボタンを用意する）

  42. VBAによる定義　①データ宣言 Private T() As Double Private U() As Double: Private V() As Double Private ID1 As Integer: Private ID2 As Integer Private KX As Integer: Private KY As Integer Private MX As Integer: Private MY As Integer Private N As Integer: Private NumLoop As Integer Private DT As Double: Private DX As Double: Private DY As Double Private R1 As Double: Private R2 As Double Private R3 As Double: Private R4 As Double Private Ndisp As Integer

  43. VBAによる定義　②初期値設定 Sub 初期値設定() With Worksheets("データ") KX = Val(.Cells(2, 1)): KY = Val(.Cells(2, 2)) DT = Val(.Cells(2, 3)): NumLoop = Val(.Cells(2, 4)) Ndisp = Val(.Cells(2, 5)) MX = KX + 1: MY = KY + 1 ReDim T(2, MX, MY), U(MX, MY), V(MX, MY) DX = 1# / KX: DY = 1# / KY R1 = 2 * DT / DX: R2 = 2 * DT / DY R3 = DT / (DX * DX): R4 = DT / (DY * DY) '　初期条件 For j = 0 To KY For i = 0 To KX T(0, i, j) = 0#: T(1, i, j) = 0# YY = DY * j U(i, j) = 50# * YY * (1# - YY) V(i, j) = 0 Next Next N = 0: ID1 = 0: ID2 = 1 End With End Sub

  44. VBAによる定義　③境界条件 Private Sub 境界条件() For j = 0 To KY T(ID1, 0, j) = 0 ‘　左側冷却 T(ID1, KX, j) = T(ID1, KX - 1, j) ‘　右側から熱流が出て行く Next For i = 0 To KX T(ID1, i, 0) = 0‘ 上・下冷却 T(ID1, i, KY) = 0 Next For i = (MX / 4 - 1) To MX / 2 - 1 T(ID1, i, 0) = 1# ‘　上部に熱源あり Next End Sub

  45. VBAによる定義　④計算実行 Private Sub 計算() 'オイラーの解法 For j = 1 To KY - 1 For i = 1 To KX - 1 T(ID2, i, j) = T(ID1, i, j) _ - R1 * U(i, j) * (T(ID1, i + 1, j) - T(ID1, i - 1, j)) _ - R2 * V(i, j) * (T(ID1, i, j + 1) - T(ID1, i, j - 1)) _ + R3 * (T(ID1, i + 1, j) - 2# * T(ID1, i, j) + T(ID1, i - 1, j)) _ + R4 * (T(ID1, i, j + 1) - 2# * T(ID1, i, j) + T(ID1, i, j - 1)) Next Next ID = ID2: ID2 = ID1: ID1 = ID End Sub

  46. VBAによる定義　⑤表示ルーチン（その１） Sub 表示() Dim CD(10) As Integer CD(0) = 5: CD(1) = 41: CD(2) = 33: CD(3) = 34: CD(4) = 36 CD(5) = 6: CD(6) = 44: CD(7) = 45: CD(8) = 46: CD(9) = 3 Worksheets("分布").Select MY = KY + 1 Application.ScreenUpdating = False Worksheets("分布").Select With Worksheets("結果") For j = 0 To KY: For i = 0 To KX .Cells(i + 1, j + 1) = T(ID1, j, i) Next: Next Umax = 1: Umin = 0 UDX = Umax - Umin

  47. VBAによる定義　⑥表示ルーチン（その２） For j = 1 To KY - 1: For i = 1 To KX - 1 If UDX = 0 Then ID = 0 Else ID = 9 * (Val(.Cells(j, i) - Umin) / UDX) End If If ID < 0 Then ID = 0 If ID > 9 Then ID = 9 Worksheets("分布").Cells(j, i).Select Selection.Interior.ColorIndex = CD(ID) Selection.Interior.Pattern = xlSolid Next: Next End With Worksheets("分布").Cells(KY + 4, KX + 4).Select Application.ScreenUpdating = True End Sub

  48. VBAによる定義　⑦ボタンのClickイベントハンドラVBAによる定義　⑦ボタンのClickイベントハンドラ Sub ボタン1_Click() 初期値設定 For i = 1 To NumLoop 境界条件 計算 If (i Mod Ndisp) = 0 Then 表示 Next End Sub

  49. 拡散の様子 ⑧結果分布図 例題による最終時刻の中央の振幅

  50. （７）その他の解法①キャビティ流れにおけるフラクショナルステップ法（７）その他の解法①キャビティ流れにおけるフラクショナルステップ法 実行条件 メッシュ数　50× 50，　レイノルド数 40, 時間刻み 0.001 時間刻み数 100，　最大収斂回数 100 図形表示の都合から，プログラムはMicrosoft C# .Netで記述 等高線：ψ，　色分け：圧力 矢印：速度ベクトル 等高線：圧力，　色分け：ψ 矢印：速度ベクトル

More Related