Tx Xx ∂ Tx () 1 1 () Xx () () ∂ t t () 1 () , t t , d dt dt t = dt dt dt d = () () () 2 Xx dx t Tx () () 2 = = = ∂ 2 c Tx ∂ = Xx () x () t c () 2 1 () t , t d d 2 dx Tx 2 dx Tx () 2 () 2

VIII.  偏微分方程式  1

2  

18．偏微分方程式と解析解  3

4  

偏微分を含む微分方程式を偏微分方程式とよぶ。多くの物理量は場の関数，すなわち時 5

間および空間の関数として表される。このため，その関係を表す方程式も時間および空 6

間の偏微分の関係として表されるので，多くの物理現象を偏微分方程式として表すこと 7

ができる。 

8

  偏微分方程式の解は境界条件によって拘束される。境界条件によって解が特殊な形で 9

書ける場合のみに解析的に解くことができる。 

10 11  

18.1  変数分離法  12

解の多変数関数が１変数関数の積で表せるときには，偏微分方程式を複数の常微分方程 13

式に分けることができる。この方法を変数分離法とよぶ。例えば１次元熱伝導方程式  14

    (18.1) 

15 16 の解を 

    (18.2) 

17

と表せるとする。これを(18.1)に代入すると  18

    (18.3) 

19

    (18.4) 

20

となる。このとき，左辺と右辺の変数を動かしてもこの式が常に成立するのは，左辺と 21

右辺がそれぞれ定数のときだけである。よって， 

22

    (18.5) 

23

    (18.6) 

24

でなくてはならない。 

25 26  

   

27

∂T x,t

( )

∂t =κ ∂²T x,t

( )

∂x²

T x,t

( )

^{=X x}

( )

^Θ

( )

^t

X x

( )

^dΘ

( )

^t

dt =κΘ

( )

t ^{d2T x}

( )

dx² 1

Θ

( )

t ^dΘ

( )

^t

dt =κ 1 X x

( )

^d

2T x

( )

dx²

1

Θ

( )

t ^dΘ

( )

^t

dt =c κ 1

X x

( )

^d

2T x

( )

dx² =c

18.2  変数変換  28

変数を別の変数に置き換えることを変数変換とよぶ。変数変換によって偏微分方程式が 29

常微分方程式に変形できることがある。熱伝導方程式において  30

    (18.11) 

31

とする。このとき， 

32

    (18.12) 

33

    (18.13) 

34

    (18.14) 

35

となるので，(18.1)に代入すると  36

    (18.15) 

37

    (18.16) 

38

となる。元の偏微分方程式に代入すると，変数x, tは消えてしまい，常微分方程式  39

    (18.17) 

40

が得られる。 

41 42   43  

19．熱伝導方程式の２つの解  44

前章の２つの方法により熱意伝導方程式の解を求める。 

45 46  

19.1  プレート冷却モデル：変数分離法による解  47

プレートモデルは熱伝導により一定の厚さのプレートの冷却を考えるモデルである。こ 48

のモデルでは両端に温度を与える，すなわち変数の値を境界条件として持つ解を求める 49

ことになる。すなわち， 

50

    (19.1) 

51

ξ= x 2 κt

∂T

∂t =∂ξ

∂t dT dξ = x

2 κ

d t

( )

^{−1 2}

dt dT

dξ =−xt^{−3 2} 4 κ

dT dξ

∂T

∂x =∂ξ

∂x dT dξ = 1

2 κt dT

dξ

∂²T

∂x² = ∂

∂x

∂T

∂x

⎛

⎝⎜

⎞

⎠⎟ = 1 2 κt

∂

∂x dT dξ

⎛

⎝⎜

⎞

⎠⎟ = 1 2 κt

∂ξ

∂x d dξ

dT dξ

⎛

⎝⎜

⎞

⎠⎟= 1 4κt

d²T dξ²

−xt^{−3 2} 4 κ

dT

dξ =κ 1 4κt

d²T dξ²

− x κt

dT dξ =d²T

dξ²

−ξdT dξ = 1

2 d²T dξ²

T

( )

0,t ^=T0

    (19.2)  52

の２つの境界値を与える。初期条件は深さ 0 以外で一定  53

    (19.3) 

54

である。プレートモデルの解は，時間が十分長く経ったときは定常解に収束することか 55

ら，(定常解＋時間依存解)という形になっているはずである。つまり， 

56

    (19.4) 

57

である。定常解は両端の温度を固定した１次元の熱伝導であることから，深さzに比例 58

59 する解 

    (19.5) 

60

である。時間依存解をT1として，これを  61

    (19.6) 

62

と書くことにする。ここで， も熱伝導方程式の解でなければならない。 に対する 63

境界条件は  64

    (19.7) 

65

    (19.8) 

66

の２つである。一方初期条件は， 

67

    (19.9) 

68

である。これは，初期の温度が深さ0以外で一定値TMとなることから来ている。また，

69

時間が無限に経った時には は0に収束する。今， が  70

    (19.10) 

71

というように深さzのみの関数と時間tのみの関数の積で表せるとする。この式を熱伝 72

導方程式(1)へ代入した式， 

73

T L,t

( )

^=TM

T z, 0

( )

^=TM

T z,t

( )

^{=T z,∞}

( )

^+T1

( )

z,t

T_∞ =T z,t

(

=∞

)

=T₀+

(

T_M −T₀

)

^z

L

T₁=

(

T_M −T₀

)

^T1′

( )

z,t

′

T₁ T₁′

T₁′

(

z=0,t

)

=0

T₁′

(

z=0,t

)

=0

T₁′

( )

z, 0 = L−z L

′

T₁ T₁′

T₁′

( )

z,t =T_z

( )

z ^Tt

( )

t

    (19.11)  74

の解は境界において0となることと，すべての点において0でない解を持つことから振 75

動的な関数でなければならないことが分かる。このため，左辺と右辺は負の一定値を持 76

つ。これを–p²で表すと， 

77

    (19.12) 

78

    (19.13) 

79

という２つの常微分方程式となる。 (19.12)の解はフーリエ級数、 

80

    (19.14) 

81

で表される。境界条件(19.7) (19.8)から，(19.14)は正弦級数でなければならない。つまり， 

82

    (19.15) 

83

である。初期条件(19.9)から、鋸歯型の波のフーリエ級数(周期 2L)となっていなくては 84

ならないことが分かる。よって， 

85

    (19.16) 

86

となる。ただし， 

87

    (19.17) 

88

である。一方，(19.13)の解は  89

   

90

である。この式は時間が十分に経つと0に収束する。また，初期条件から， 

91

    (19.18) 

92

でなければならない。よって求める温度は，フーリエ級数を用いて  93

    (19.19) 

94

と求められる。有限の時間ではLを十分大きく取れば，半無限体冷却モデルの解も表す 95

ことができる。下図は，(19.19)によるプレート冷却モデルの解である。0 Maのときの解 96

には，不連続関数をフーリエ級数で表したことによるギプスの現象が見られる。 

97

1 κ

1 T_t

∂Tt

∂t = 1 T_z

∂²T_z

∂z²

d²T_z

dz² =−p²T_z 1

κ dT_t

dt =−p²T_t

T_z =

(

A_ncosp_nz+B_nsinp_nz

)

n=0

∑

∞

A_n =0

 T_z = 2

nπsin nπz L

⎡

⎣⎢

⎤

⎦⎥

n=1

∑

∞

p_n = nπ L

T_t =C_nexp⎡⎣−p_n²κt⎤⎦

C_n =1

 

T =T₀ +

(

T_M −T₀

)

^z

L + 2

nπsin nπz L

⎡

⎣⎢

⎤

⎦⎥exp −n²π²κt L²

⎡

⎣⎢ ⎤

⎦⎥

n=1

∑

∞

⎡

⎣⎢ ⎤

⎦⎥

98  

99  

図 19.1  フーリエ級数によるプレート冷却モデルの解 

100 101  

19.2  半無限体冷却モデル：変数変換法による解  102

  半無限体冷却モデル  (Half-Space Cooling model: HSCM)は，プレートを無限に下方 103

へ続くマントルが上面から冷却された部分であると考えるモデルである。境界条件は  104

    (19.21) 

105

    (19.22) 

106

であり，この条件はフーリエ級数では表すことができない。ここでは，変数変換を利用 107

108 する。 

    (19.23) 

109

と置くと，(18.17)の解は  110

    (19.24) 

111

T z

(

=0,t

)

^=T0

T z

(

=∞,t

)

^=TM

φ=dT dξ

φ=C₁e^−ξ2

と表される。(19.23)より，この式をxで積分すると  112

    (19.25) 

113

変数変換した境界条件  114

    (19.26) 

115

    (19.27) 

116 117 より， 

    (19.28) 

118

    (19.29) 

119

となる。よって， 

120

    (19.30) 

121

    (19.31) 

122

となる。すなわち，温度は  123

    (19.32) 

124

と表される。右辺の積分の部分  125

    (19.33) 

126

は誤差関数とよばれる。誤差関数はガウス関数を積分した関数で，0から1までの値を 127

とる。誤差関数はxが無限大のときを除いて解析的に解けないので，数値積分により求 128

める必要がある。 

129 130  

T

( )

ξ ^=C0+C₁ e^{− ′ξ2}dξ′

0

∫

ξ

T

(

ξ=0

)

^=T

( )

^{0,t =T}0

T

(

ξ=∞

)

^=T

( )

^{∞,t =T}M

T₀=C₀+C₁⋅0

T_M =C₀+C₁ e^{− ′ξ2}dξ′

0

∫

∞ ^=T0+C1 ₂^π

C₀=T₀

C₁=

(

T_M −T₀

)

²_π

T

( )

ξ ^=T0+

(

^{TM −T0}

)

²_π

∫

₀^{ξe− ′ξ2dξ′}

erf

[ ]

ξ ^{= 2}_π

∫

₀^{ξe− ′ξ2dξ′}

131  

図 19.2  HSCM による地殻熱流量とプレートの温度の時間変化：地殻熱流量の観測値は

132

Davies (2013)によるグリッドデータから作成 

133 134   135  

20．酔歩運動と拡散  136

  拡散とは，酔歩運動(ランダムウォーク)によって不均質が物質中に広がっていく過程 137

である。熱伝導も原子の不規則な衝突により原子振動が拡散していく現象であり，拡散 138

方程式と熱伝導方程式は同じ形で表すことができる。 

139 140  

20.1  酔歩運動とフォッカー・プランク方程式  141

  粒子がdtの間にdx だけ進むとする。ここでdx はランダムに決まるとする。このよう 142

な場合，t における粒子の位置を決めるにはt – dt における粒子の位置情報だけが必要 143

であり，それより前の情報は必要がない。このような過程をマルコフ過程とよぶ。粒子 144

の密度f(x, t)を求めることを考える。ここで，粒子の密度は単に長さ辺りの粒子数とす

145

る。時刻t–d_tにおける粒子の密度がf(x, t – dt)で与えられているとき，x – dxからxへ粒 146

子が動く確率をP(x – dx, x)とする。そのとき，f(x, t)は  147

    (20.1)  148

と表すことができる。ここでfをx, tのまわりでテイラー展開したもの  149

     

150

      (20.2) 

151

を右辺に代入すると， 

152

     

153

            (20.3) 

154

         

155

         

156

         

157

              (20.4) 

158

         

159

         

160

         

161

          (20.5) 

162

となる。ここで，P(x – dx, x)は確率なので  163

    (20.6) 

164

f x,t

( )

⁼

∫

_−∞^{∞ f x}

(

^{−δx,t−δt}

)

^{P x}

(

^−δx,x

)

^d

( )

^δx

f x

(

−δx,t−δt

)

^{= f x,t}

( )

^{− ∂}_∂x^fδx+1

2

∂²f

∂x²δx²−!

− ∂f

∂t δt+1 2

∂² f

∂t² δt²−!+ ∂²f

∂t∂xδtδx−!

f x,t

( )

^{= f x,t}

( )

^{− ∂}_∂x^fδx+1 2

∂²f

∂x²δx²−!

⎡

⎣⎢

−∞

∫

∞ ^{− ∂}_∂t^fδt+1₂^∂_∂t^{22f δt2−!}

+ ∂²f

∂t∂xδtδx−!⎤

⎦⎥⋅P x

(

−δx,x

)

^d

( )

^δx

=

∫

_−∞^∞ f x,t

( )

^{P x}

(

^−δx,x

)

^d

( )

^δx

− ∂f

∂xδxP x

(

−δx,x

)

^d

( )

^δx

−∞

∫

∞ ⁺ −∞¹₂^∂_∂x^{2f2δx2P x}

⁽

^−δx,x

⁾

^d

^{( )}

^δx

∫

∞

− ∂f

∂t δtP x

(

−δx,x

)

^d

( )

^δx

−∞

∫

∞ ⁺ −∞¹₂^∂_∂t^{22fδt2P x}

⁽

^−δx,x

⁾

^d

^{( )}

^δx

∫

∞

+ ∂² f

∂t∂xδtδxP x

(

−δx,x

)

^d

( )

^δx

−∞

∫

∞ ^−!

= f x,t

( ) ∫

_−∞^{∞ P x}

(

^−δx,x

)

^d

( )

^δx

− ∂f

∂x

∫

_−∞^∞δxP x

(

−δx,x

)

^d

( )

^{δx +1}₂^∂_∂x^2f2 −∞^{δx2P x}

⁽

^−δx,x

⁾

^d

^{( )}

^δx

∫

∞

− ∂f

∂t δt

∫

_−∞^∞ P x

(

−δx,x

)

^d

( )

^{δx +1}₂^∂_∂t^{22f δt2} −∞^{P x}

⁽

^−δx,x

⁾

^d

^{( )}

^δx

∫

∞

+ ∂² f

∂t∂xδt

∫

_−∞^∞δxP x

(

−δx,x

)

^d

( )

^{δx −!}

P x

(

−δx,x

)

^d

( )

^δx

−∞

∫

∞ ⁼¹

    (20.7)  165

    (20.8) 

166

である。それぞれ，確率の合計，平均値，分散を表す。このことに注意すると，(20.5) 167

168 は， 

   

169

            (20.9) 

170

となる。ここで高次の項を無視して，さらに移項すると  171

    (20.10) 

172

となり， 

173

    (20.11) 

174

が得られる。ここでもd_t→0とすると，<d_x>も0となるので， 

175

    (20.12) 

176

が得られる。この式はフォッカー・プランク方程式(Fokker-Plank  equation)とよばれ 177

る。ここで， 

178

    (20.13) 

179

は流速， 

180

    (20.14) 

181

は拡散係数を表す。つまり，この式は移流拡散方程式  182

    (20.15) 

183

を表す。左右へ動く確率が等しい，つまり， 

184

    (20.16) 

185

ときには  186

δxP x

(

−δx,x

)

^d

( )

^δx

−∞

∫

∞ ^{= δx}

δx²P x

(

−δx,x

)

^d

( )

^δx

−∞

∫

∞ ^{= δx2}

f x,t

( )

^{= f x,t}

( )

⁻ δx ∂f

∂x + δx² 2

∂² f

∂x² − ∂f

∂t δt+1 2

∂²f

∂t² δt² + δx ∂²f

∂t∂xδt+O

(

δx³,δt³

)

∂f

∂t δt=− δx ∂f

∂x + δx² 2

∂² f

∂x² +1 2

∂² f

∂t² δt²+ δx ∂²f

∂t∂xδt

∂f

∂t =− δx δt

∂f

∂x + δx² 2δt

∂² f

∂x² +1 2

∂² f

∂t² δt+ δx ∂²f

∂t∂x

∂f

∂t =− δx δt

∂f

∂x + δx² 2δt

∂² f

∂x²

v= δx δt

D= δx² 2δt

∂f

∂t =−v∂f

∂x +D∂² f

∂x² δx =0

    (20.17)  187

となり，拡散方程式を表す。   

188 189  

20.2  酔歩運動と拡散方程式の解  190

  ここでは，初期条件が原点に拡散する不均質が集中している場合を考える。この場合，

191

初期濃度はディラックのデルタ関数により， 

192

    (20.18) 

193

と表される。拡散方程式の解は，(19.12) (19.13)より，フーリエ級数をフーリエ変換に置 194

き換えた， 

195

    (20.19) 

196

である。ここで，t = 0でfがデルタ関数となるための条件は， 

197

   

198

である。よって， 

199

    (20.20) 

200

となる。ここで，17 章の(17.31)より， 

201

    (20.21) 

202

を用いると， 

203

    (20.22) 

204

    (20.23) 

205

と置き換えれば良いので，求める解は， 

206

    (20.24) 

207

となる。 

208

  数値シミュレーションを行って，酔歩運動により拡散現象が再現できるか試した結果 209

が図 20.1 である。粒子は最初に原点x = 0にN0個密集しているとする。粒子の番号を 210

∂f

∂t =D∂² f

∂x²

f x,0

( )

^=δ

( )

^x

f x,t

( )

⁼_π¹

_∫

₀^{∞A p}

( )

^cospxexp_⎡⎣^−p2Dt_⎤⎦^dp

A p

( )

⁼¹

f x,t

( )

⁼_π¹

_∫

₀^∞cospxexp_⎡⎣^−Dtp2_⎤⎦^dp

exp⎡⎣−ax²⎤⎦cosbx dx

0

∫

∞ ^{= 1}₂ ^π_a^exp⎡_⎣^⎢−_4a^{b2 ⎤}_⎦^⎥

a= Dt b=x

f x,t

( )

^{= 1}_{2 π}¹

Dtexp − x² 4Dt

⎡

⎣⎢ ⎤

⎦⎥

iとすると，初期条件で  211

    (20.25) 

212

である。ここで<dx²> = 1となるような正規乱数を発生させ， 

213

    (20.26) 

214

を計算する。ただしdt = 1とする。このとき，粒子の座標を区間Dxに区切って区間ごと 215

の粒子数をカウントする。ところで，粒子数と拡散係数はそれぞれ， 

216

    (20.27) 

217

    (20.28) 

218

であるので，その解は  219

    (20.29) 

220

となる。ただし，nはタイムステップを表す。この式は平均０分散nの正規分布を表す。 

221 222  

223  

図 20.1  ランダムウォーク粒子計算と拡散方程式の解の比較 

224 225  

x_i=0

x_iⁿ⁺¹=x_iⁿ+δx_i

f x,t

( )

−∞

∫

∞ ^dx=N0

D= 1 2

f x,

(

tⁿ

)

⁼ ₂^N_π^o_n^exp"_#^$−x_2n^2%_&^'