Excelで定積分・微分方程式

８．４偏微分方程式の数値解法
［偏微分の復習（１）］
f ( x, y)  x2  xy  y 2 のとき
df ( x, y)
dy 
dy

 2x   y  x   2 y
dx
dx 
dx

この x と y を xy 座標系の座標値とすると，
お互いに直交する値であるから， x と y は独立。
dy dx は意味がない。
そこで
お互いに相手の変数を定数とみなして微分する。
偏微分と呼ぶ。
［記号］ d の替わりに  を用いる。
f ( x , y )
x
f ( x, y)  x  xy  y
2
2
偏微分の例
y を定数とみなして
f ( x, y )
 2 x  y　x
x を定数とみなして
f ( x, y )
 x  2y
x
 を見慣れないからといって恐れないこと。常微分より簡単だよ。
［偏微分の復習（２）］
［常微分と偏微分の関係］
df ( x, y ) f ( x, y ) f ( x, y ) dy


dx
x
y dx
［媒介変数表示］
x(u, v), y(u, v)が，ともに u, v に関して偏微分できるとき，
f f x f y


,
u x u y u
f f x f y


v x v y v
以上の関係は，いたるところで出てくるので重要！！
［ベクトル解析の復習①］
［勾配］
grad( f ) 
f
f
f
i
j k
x
y
z
i , j , k は，それぞれ x, y , z 方向の単位ベクトル
i



 j k
という演算子（ナブラ）を導入すると
x
y
z
f  grad( f )
２次元の地形図を思い浮かべるとよい。
イメージをつかむために


j
x
y
スカラー値 f を標高値とすると， f は傾斜方向ベクトルである。
２次元の▽（ナブラ）は   i
f
f  i
f
f
j
x
y
絶対値：
y
f
x   f x 
2
 f y 

 f x 
方向：   Tan1 

x
f y
f x
2
［ベクトル解析の復習③］


［発散］ベクトル場 A  Ax , Ay , Az に対して
div A    A 
Ay
Ax
A

 z
x
y
z
ベクトル場 A がスカラー f の勾配のとき
簡単のために２次元の数値地図を考え，
左のメッシュでは Ax  0, Ay 不変とする
と， div A が正のとき，傾斜方向は広が
り，負のときは狭まる形になる。
2 f 2 f 2 f
div( grad f )    f  2  2  2
x
y
y
ここで，   f  2 f  f とあらわし
これを f のラプラシアンという。
Ax
0
x
Ax
0
x
Ax
0
x
（正のとき凹型地形，負のとき凸型地形）
（あるいは，隣り合うベクトルの広がり具合）
［ベクトル解析の復習④］


［回転］ベクトル場 A  Ax , Ay , Az に対して
rot A    A
 Az Ay   Ax Az   Ay Ax 
i  
k
 



 j  
z   z
x   x
y 
 y
なお，
i

rot A    A  
 x
 Ax
j

y
Ay
k

 と表記されることもある。
z 
Az 
［ベクトル解析の復習⑤］
［公式］
公式を丸暗記しないこと。意味を理解する。できなければ，どの資料に書い
てあるかを記憶し，使うときに参照する。使っているうちに慣れてくる。
① grad( f  g )  ( f  g )  f  g
 
 


② div( A  B)    ( A  B)    A    B
 
 


③ rot( A  B)    ( A  B)    A    B
④ div grad f    f   2 f

⑤ rot grad f    f  0
  f    f 
(r o gt r a fd) x        0 となり，
y  z  z  y 
その他の成分も同様に 0 となる。
［ベクトル解析の復習⑥］


⑥ div rot A    (  A)  0
   Az Ay    Ax Az    Ay Ax 
  

div rot A  



  
x  y
z  y  z
x  z  x
y 
2
2
 2 Az  Ay  2 Ax  2 Az  Ay  2 Ax






0
xy xz yz yx zx zy




2
⑦ rot rot A    (  A)  (  A)   A

  Ay Ax    Ax Az 
  
(rot rot A) x  



y  x
y  z  z
x 
  Ax Ay Az    2 Ax  2 Ax Ax
  
 



 2
2
2
x  x
y
z   x
y
z


2
 ((  A)) x  ( A) x



［ベクトル解析の復習⑦］
⑧ grad( fg )  ( fg )  gf  f g




⑨ div( fA)    ( fA)  (f )  A  f   A




⑩ rot( fA)    ( fA)  (f )  A  f   A
 
 
  

⑪ div( A  B)    ( A  B)  (  A)  B  A  (  B)
 
 



 


⑫ rot( A  B)    ( A  B)  ( B  ) A  ( A  ) B  A (  B)  B(   A)
 
 




 

⑬ grad( A  B)  ( A  B)  ( B  ) A  ( A  ) B  B  (  A)  A  (   B)
８．４．１数値解法における差分の方法
（１）差分の種類
常微分の場合と同様，以下のような方法がある。
①前進差分（forward difference）
②中心差分（centered difference）
③後退差分（backward difference）
（２）前進差分
微小格子を考え，以下のように近似する。
u ui 1, j  ui , j

x
h
u ui , j 1  ui , j

y
h
 2u ui  2, j  2ui 1, j  ui , j

2
x
h2
 2u ui 1, j 1  ui 1, j  ui , j 1  ui , j

xy
h2
 2u ui , j  2  2ui , j 1  ui , j

2
y
h2
ui , j  2
ui , j 1
ui 2, j ui 1, j
ui , j
ui , j 1
ui , j 2
h
ui 1, j ui  2, j
h
（３）後退差分
微小格子を考え，以下のように近似する。
u ui , j  ui 1, j

x
h
u ui , j  ui , j 1

y
h
 2u ui , j  2ui 1, j  ui  2, j

2
x
h2
 2u ui , j  ui 1, j  ui , j 1  ui 1, j 1

xy
h2
 2u ui , j  2ui , j 1  ui , j  2

2
y
h2
ui , j  2
ui , j 1
ui 2, j ui 1, j
ui , j
ui , j 1
ui , j 2
h
ui 1, j ui  2, j
h
（４）中心差分
微小格子を考え，以下のように近似する。
u ui 1, j  ui 1, j

x
2h
u ui , j 1  ui , j 1

y
2h
 2u ui 1, j  2ui , j  ui 1, j

2
x
h2
 2u ui 1, j 1  ui 1, j 1  ui 1, j 1  ui 1, j 1

xy
4h 2
 2u ui , j 1  2ui , j  ui , j 1

2
y
h2
ui , j  2
ui , j 1
ui 2, j ui 1, j
ui , j
ui 1, j ui  2, j
ui , j 1
ui , j 2
h
中心差分は，前進差分から後退差分を差し引いて，１／２にした形である。
h
（５）差分式使用上の注意点
どの方法を選択するかは，
①与えられた条件，
②問題の性格，
③モデル化のしやすさ
等に応じて決めること。
流体力学の分野では，空間成分に関して，速度成
分が正であれば後退差分を，負であれば前進差
分を使用する風上差分も用いられている。
８．４．２物理現象と偏微分方程式
（１）考え方
①偏微分方程式は，物理現象を支配する法則等により
体積無限小の極限で導出される。
②有限な体積要素における法則により，基礎式が直接
導出されることも多い。
（例：コントロールボリューム法，直接差分法）
２階偏微分方程式の一般式
（２）２階偏微分方程式の分類

 2u
 2u
 2u
u u 
A 2 B
 C 2  F  x, y, u, , 
x
xy
y
x y 

［分類］
① B 2  4 AC  0 のとき楕円型
［例］  2u
 2u
 2 0
2
x
y
（ラプラス方程式：Laplace’s equation）
 2u  2u
 2  f x, y  （ポアソン方程式：Poisson’s equation）
2
x
y
② B 2  4 AC  0 のとき放物型
2

u

u
［例］

 f  x, t 
t
x 2
2
③ B  4 AC  0 のとき双曲型
（拡散方程式など）
2
2

u

u
2
［例］

c
 f x, t  （波動伝播など）
2
2
t
x
 2u  2u
 2 0
2
x
y
（ラプラス方程式）
 2u  2u
 2  f x, y  （ポアソン方程式）
2
x
y
［シートの定義］
（３）楕円型のプログラム例
VBAによる定義
①データ宣言
Private
Private
Private
Private
Private
Private
Private
Private
Private
Private
Private
Private
ポアソン As Boolean
KX As Integer
KY As Integer
NN As Integer
EPS As Double
U() As Double
ID1 As Integer
ID2 As Integer
F() As Double
UMax As Double
UMin As Double
UDX As Double
②初期条件の設定と境界条件設定
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
③計算本体（その１）
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
④計算本体（その２）
' 収斂計算
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
⑤計算本体（その３）
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
⑥データの設定
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
⑦結果の設定
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
⑧ボタンのClickイベントハンドラ
Sub ボタン1_Click()
‘ラブラスの方程式のためのClickイベントハンドラ
データ設定
ポアソン = False
計算 KX, KY, NN, EPS
結果設定
End Sub
Sub ボタン3_Click()
‘ポアソンの方程式のためのClickイベントハンドラ
データ設定
ポアソン = True
計算 KX, KY, NN, EPS
結果設定
End Sub
⑨結果を分布図としてセルの色で表示する処理
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
⑩結果分布図
ラプラスの方程式の結果
(Umax=15.35, Umin=0.00)
ポアソンの方程式の結果
(Umax=14.93, Umin=0.03)
ラプラスの方程式の結果の方が分布として
滑らかであることに注意されたい。
電磁場における支配方程式
（マクスウェルの方程式）
（４）電磁場の偏微分方程式
（ガウスの法則）
ポアソンの方程式
D  
B
 E 
0
（電磁誘導の法則）
t
B  0
（磁束密度湧き出しなしの法則）ラプラスの方程式
D
 H 
 J （アンペールの法則）
t
：電束密度
D
E
H
B
J

：電場の強さ
：磁場の強さ
：磁束密度
：電流密度
：電荷
基本式が同じでも，
条件によって
偏微分方程式が異なることを
ここでは
確認する。
（４）電磁場の偏微分方程式
誘電率  ，透磁率  ，電荷 
を用いた電束密度と電場強さ，磁束密度と
磁場の強さ，電流密度と電場の強さの関係
D  E , B  H , J  E
これらをマックスウェルの方程式に代入して
 E  
H
E
 0,   H  
 E
t
t
右の式を時間で偏微分して，左の式を代入すると
H
2 E
E

 2 
t
t
t
2
2 E
E
1
1

E
 2 
      E        E    2 E 
t
t





条件で変わる偏微分方程式
c
1

同じ方程式でも，
条件によって偏微分方程式が異なる
として整理すると
 2 H  E
2 2


c
 E
2
t
 t
真空中では
 0
2 H
2 2

c
 E
2
t
（双曲型：真空中の電場伝播）
電場の強さの
時間変化が小さいとき
 E
 c 2 2 E
 t
（放物線型：拡散方程式）
定常状態のとき元の式
H
 E  
 0   E  0
t
電位  を用いて
E  
 2  


（楕円型：ポアソン方程式）
（５）波動方程式
１次元（平面波）
2
 2u
2  u

2
t
x 2
２次元（円柱波）
2
2
 2u
2  u
2  u
 x 2 y 2
2
t
x
y
（梁の振動：１次元）
4
 2u
2  T

2
t
x 4
３次元（球面波）
2
2
2
 2u
2  u
2  u
2  u
 x 2  y 2 z 2
2
t
x
y
z
極座標による２次元波動方程式（太鼓など円形の膜の振動のとき）
2
2

 2u

u
1

u
1

u
2

 c  2 

2
2 
t
r r r  
 r
１次元波動方程式の解釈
１次元（平面波）
2
 2u
2
2  u
u j 1  2u j  u j 1 


2
2
2
t
x
(x)
隣り合う両節点の変位差を変位とみなせば，
加速度が変位に比例するバネと同じ
u
j 1
u j 1  u j
 u j   u j  u j 1 
 u j 1  2u j  u j 1
u j  u j 1
x
x
u
 v,
t
プログラム例
2
v
2  u

t
x 2
とおいて，以下の差分方程式としてプログラミングする。
ui 1  ui  vt
vi 1  vi  t 
2
(x)
2
u
j 1
 2u j  u j 1 
［シートの定義］（以下のように３種類のシートとコマンドボタンを用意する）
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
‘ 端点での誤差を少なくするために
‘ このように設定している
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
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
中央の振幅
⑨結果グラフ
両端での速度・位置を固定しているので，
振幅が次第に小さくなる。
T=100までの弦の変化
αの値を変化させて，
色々確かめてみよう。
また，αが大きすぎると
振動状態になってしまうことも
確認すること。
（６）保存量の偏微分方程式
熱拡散方程式
境界面を移動する流束ベクトル
 T  T  T 
T
   2  2  2   Q
t
y
z 
 x
2
2
J  J x , J y , J z 
2
 J z 
Jz  
z
 z 
Q  0 の条件による定常状態では
 J y
J y  
 y
 2T  2T  2T
 2  2 0
2
x
y
z
Jx
z
 J 
J x   x x
 x 
Jy
y
x
Jz

 y

プログラム例
［シートの定義］（以下のように３種類のシートとコマンドボタンを用意する）
VBAによる定義
①データ宣言
Private
Private
Private
Private
Private
Private
Private
Private
Private
Private
Private
T() As Double
U() As Double: Private V() As Double
ID1 As Integer: Private ID2 As Integer
KX As Integer: Private KY As Integer
MX As Integer: Private MY As Integer
N As Integer: Private NumLoop As Integer
DT As Double:
DX As Double: Private DY As Double
R1 As Double: Private R2 As Double
R3 As Double: Private R4 As Double
Ndisp As Integer
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
VBAによる定義
②初期値設定
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
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,
- R2 * V(i, j) * (T(ID1, i, j +
+ R3 * (T(ID1, i + 1, j) - 2# *
_
+ R4 * (T(ID1, i, j + 1) - 2# *
Next
Next
ID = ID2: ID2 = ID1: ID1 = ID
End Sub
VBAによる定義
④計算実行
j) - T(ID1, i - 1, j)) _
1) - T(ID1, i, j - 1)) _
T(ID1, i, j) + T(ID1, i - 1, j))
T(ID1, i, j) + T(ID1, i, j - 1))
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
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
VBAによる定義
⑦ボタンのClickイベントハンドラ
Sub ボタン1_Click()
初期値設定
For i = 1 To NumLoop
境界条件
計算
If (i Mod Ndisp) = 0 Then
Next
End Sub
表示
拡散の様子
⑧結果分布図
例題による最終時刻の中央の振幅
（７）その他の解法
①キャビティ流れにおけるフラクショナルステップ法
図形表示の都合から，プログラムはMicrosoft C# .Netで記述
実行条件メッシュ数 50× 50，レイノルド数 40, 時間刻み 0.001
時間刻み数 100，最大収斂回数 100
等高線：ψ，色分け：圧力
矢印：速度ベクトル
等高線：圧力，色分け：ψ
矢印：速度ベクトル
（７）その他の解法
②キャビティ流れにおける流れ関数・過度法
図形表示の都合から，プログラムはMicrosoft C# .Netで記述
実行条件メッシュ数 50× 50，レイノルド数 40, 時間刻み 0.001
時間刻み数 100，最大収斂回数 100，加速度パラメータ 1.0
等高線：ψ，色分け：ω
等高線：ω，色分け：ψ
図形表示の都合から，プログラムは
Microsoft C# .Netで記述
（等高線，色分けともに：ψ）
（７）その他の解法
③円柱回りの流れ

Download Report