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

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

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

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

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

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

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

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

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

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

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

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

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

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

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 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 End With End Sub

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

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 イベントハンドラ 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 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 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 イベントハンドラ Sub ボタン 1_Click() 初期値設定 For i = 1 To NumLoop 境界条件 計算 If (i Mod Ndisp) = 0 Then 表示 Next End Sub

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

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

51 （７）その他の解法 ② キャビティ流れにおける流れ関数・過度法 図形表示の都合から，プログラムは Microsoft C#.Net で記述 等高線： ψ ， 色分け： ω 等高線： ω ， 色分け： ψ 実行条件 メッシュ数 50× 50 ， レイノルド数 40, 時間刻み 0.001 時間刻み数 100 ， 最大収斂回数 100 ， 加速度パラメータ 1.0

52 （７）その他の解法 ③円柱回りの流れ 図形表示の都合から，プログラムは Microsoft C#.Net で記述 （等高線，色分けともに： ψ ）