拡散方程式の演習 - Moriturus te saluto

PDE勉強し始めたので一番簡単なやつからやっつける。

とりあえず鼻息荒げにカテゴリ2種類も作ったので、この記事だけでおしまいみたいな事にはしたくないですね。

熱伝導方程式のモデル

例題

長さ $L = 2$ [m]の銅線を温度（ $T_0 = 10^\circ \mathrm{C}$ ）の環境中に置く。銅線は長い間置かれており、温度変化のない環境と変わらない温度であると仮定する。

いま、時間 $t = 0$ においてこの銅線の左端を温度 $T_1 = 0^\circ \mathrm{C}$ , 右端を温度 $T_2 = 50^\circ \mathrm{C}$ の熱源でつないだ。

この時の銅線の場所ごとに時間温度変化を調べたい。

ヒント：この時間による温度の変化は放物型偏微分方程式で記述できる。

以上に基づいて3つの条件、

を決める。

熱伝導方程式

一次元空間 $x$ における熱伝導を考えたとき、以下の様な関係式が成り立つ(PDE)。

$u_t = \alpha^{2} u_{xx} \hspace{20.0em} 0 < x < L \hspace{20.0em} 0 < t < \infty$

$u_t$ は時間による温度の変化率、 $u_{xx}$ は $u(x, t)$ の凹面である。

$u_{xx} \approx \frac{1}{\Delta x^{2}} (u(x + \Delta x, t) - 2u(x, t) + u(x - \Delta, t))$

これを次のように書き換えると、隣接する地点の平均と $u(x, t)$ の差分を取っていると解釈できる。

$u_{xx} \approx \frac{2}{\Delta x^{2}} \biggl( \frac{u(x + \Delta x, t) + u(x - \Delta x, t)}{2} - u(x, t) \biggr)$

境界条件

問題文の条件を当て込むだけ。

$\begin{eqnarray} Boundary\hspace{4.0em}Conditions=\left\{\begin{array} u(0, t) = T_1 \\ u(L, t) = T_2 \end{array}\right. \hspace{6.0em} 0 < t < \infty \end{eqnarray}$

初期条件

$t = 0$ において $u(x, 0)$ はあまねく温度 $T_0 = 10^\circ \mathrm{C}$ である。

数値計算

これらの条件から $u(x, t)$ を求めるのが目的であるが、ここまでの情報ですでに数値計算によって拡散方程式の定性的な動態を観察することはできるので、さしあたりそちらで様子を見てみる。ちなみにコードはPython+matplotlib+NumPy

import numpy as np
from matplotlib import pyplot as plt
from matplotlib import animation as animation

# Calculate parabolic PBE
def heat_diffusion(xs, ys, dx):
    for i, x in enumerate(xs):
        if i == 0 or i == len(xs) - 1:
            continue
        ys[i] += ((ys[i-1] + ys[i+1])/2. - ys[i])
    return ys

def main():
    # configuration
    dx = 1. / 200.
    max_time = 10000

    # initialization
    xs = np.arange(0., 2., dx)
    ys = np.copy(xs)
    ys.fill(10)

    # boundary condition
    ys[0] = 0.
    ys[len(ys)-1] = 50.

    for i in xrange(0, max_time):
        ys = heat_diffusion(xs, ys, dx)
        if i % (max_time/10) == 0:
            plt.plot(xs, ys, label="Time="+str(i))

    plt.legend(loc="upper left")
    plt.xlabel("position [m]")
    plt.ylabel("temperature [C]")
    plt.title("Heat Diffusion")
    plt.show()

if __name__ == "__main__":
    main()

f:id:xef:20130127130457g:plain

Reference

Partial Differential Equations for Scientists and Engineers, Stanley J. Farlow