1次元熱伝導方程式をPythonで解く

はじめに

熱伝導方程式は、物質内の熱の伝わり方を記述する基本的な偏微分方程式の一つです。本記事では、一次元の熱伝導方程式をPythonを用いて数値的に解く方法と、その解析解を求める方法を紹介します。さらに、得られた数値解と解析解をアニメーションとして可視化します。

熱伝導方程式とは？

一次元の熱伝導方程式は次のように表されます：

$$ \frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2} $$

ここで、$u(x,t)$ は時刻 $t$ における位置 $x$ での温度、$\alpha$ は熱拡散率です。

熱伝導方程式がなぜその形の式なのか？と思った方は、下記にざっくりとした説明を載せたのでそちらを参考にしてください。フォッカー-プランク方程式で、ドリフト項がゼロで、拡散係数が一定である、と仮定したのが熱伝導方程式です。

タイムスケールの計算方法

一次元の熱伝導方程式において、空間と時間のオーダーを近似するためには、無次元化を通じてオーダー解析を行います。

一次元熱伝導方程式

まず、一次元の熱伝導方程式を確認します。

空間スケール $ L $ と時間スケール $ T $ を用いて、変数を無次元化します。

x' = \frac{x}{L}, \quad t' = \frac{t}{T}, \quad u' = \frac{u}{U}

この変換を用いると、元の方程式は次のように書き換えられます。

\frac{U}{T} \frac{\partial u'}{\partial t'} = \alpha \frac{U}{L^2} \frac{\partial^2 u'}{\partial x'^2}

次に、すべての項を無次元化します。

オーダー解析

各項の次元解析を行い、無次元数を導出します。

\frac{U}{T} \sim \frac{\alpha U}{L^2}

これを整理すると、

T \sim \frac{L^2}{\alpha}

となります。ここから、時間スケール $ T $ と空間スケール $ L $ の関係が導かれます。つまり、

空間スケール $ L $ が決まれば、時間スケール $ T $ は $ T \sim \frac{L^2}{\alpha} $ で近似できます。

このオーダー解析により、熱伝導の問題において空間スケールと時間スケールのオーダーを抑えてから、具体的な計算をしましょう。

初期条件と境界条件

今回は、次のような初期条件と境界条件を設定します：

初期条件：方形波（最大値5）
境界条件：固定境界条件（両端は常に0度）

具体的には、次のように設定します：

u(x, 0) = 
\begin{cases} 
   \frac{5}{L/2} x & \text{if } 0 \leq x \leq \frac{L}{2} \\
   \frac{5}{L/2} (L - x) & \text{if } \frac{L}{2} < x \leq L 
\end{cases}

ここで、$L$ は棒の長さ（今回は10）です。

解析解の導出

初期条件が与えられた場合の解析解は、フーリエ級数展開を用いて次のように求めることができます：

1) 初期条件のフーリエ係数の計算

まず、初期条件をフーリエ級数展開します：

u(x, 0) = \sum_{n=1}^{\infty} b_n \sin\left(\frac{n \pi x}{L}\right)

離散フーリエ級数展開なので、偶関数であればコサインが残り、今回のように $x$ について奇関数であれば、サインの項だけが残ります。無限個の足し算はできませんので、数値的に解くためには、有限個の足し算で止めることになります。

ここで、フーリエ係数 $b_n$ は次のように計算されます：

b_n = \frac{2}{L} \int_0^L u(x, 0) \sin\left(\frac{n \pi x}{L}\right) \, dx

この積分を行うと、次の結果が得られます：

b_n = \frac{20}{(n \pi)^2} (1 - (-1)^n)

2) 温度分布の時間発展

得られたフーリエ係数を用いて、温度分布の時間発展を表すフーリエ級数を構築します：

$$ u(x, t) = \sum_{n=1}^{\infty} b_n \sin\left(\frac{n \pi x}{L}\right) \exp\left(-\alpha \left(\frac{n \pi}{L}\right)^2 t\right) $$

これにより、任意の時刻 $t$ における温度分布 $u(x, t)$ を計算できます。

ここで、t=0での初期条件から指数関数的な時間発展はどう導出したのか？と疑問に思った人は次の各モードの時間発展も読んでください。

各モードの時間発展

フーリエ級数展開を用いると、温度分布 $u(x, t)$ を次の形に書き直せます：

u(x, t) = \sum_{n=1}^{\infty} u_n(t) \sin\left(\frac{n \pi x}{L}\right)

ここで、$u_n(t)$ は時間依存の係数です。これを熱伝導方程式に代入します。

\frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2}

すると、

\sum_{n=1}^{\infty} \frac{d u_n(t)}{d t} \sin\left(\frac{n \pi x}{L}\right) = - \alpha \sum_{n=1}^{\infty} u_n(t) \left(\frac{n \pi}{L}\right)^2 \sin\left(\frac{n \pi x}{L}\right)

同じフーリエモードについて比較すると、

\frac{d u_n(t)}{d t} = -\alpha \left(\frac{n \pi}{L}\right)^2 u_n(t)

微分方程式の解

この微分方程式は、各フーリエモードに対して次の形をしています：

\frac{d u_n(t)}{d t} = -\alpha \left(\frac{n \pi}{L}\right)^2 u_n(t)

これは線形の常微分方程式で、解は次のように指数関数的な減衰を示します：

u_n(t) = u_n(0) \exp\left(-\alpha \left(\frac{n \pi}{L}\right)^2 t\right)

ここで、初期条件 $t = 0$ における $u_n(0)$ は、初期温度分布のフーリエ係数 $b_n$ です。したがって、

u_n(t) = b_n \exp\left(-\alpha \left(\frac{n \pi}{L}\right)^2 t\right)

温度分布の時間発展

これにより、任意の時刻 $t$ における温度分布 $u(x, t)$ は次のように表されます：

u(x, t) = \sum_{n=1}^{\infty} b_n \sin\left(\frac{n \pi x}{L}\right) \exp\left(-\alpha \left(\frac{n \pi}{L}\right)^2 t\right)

これが、初期条件をフーリエ級数展開で求めた後の、温度分布の時間発展の一般的な解法です。

有限差分法（Finite Difference Method, FDM）を用いた数値的解法

熱伝導方程式を有限差分法（Finite Difference Method, FDM）を用いて離散化して解きます。
有限差分法とは、偏微分方程式を差分方程式に変換して、数値的に解くための方法です。

空間微分の離散化

空間微分 $\frac{\partial^2 u}{\partial x^2}$ を中央差分法で離散化します。これにより、次のようになります：

\frac{\partial^2 u}{\partial x^2} \approx \frac{u(x+\Delta x, t) - 2u(x, t) + u(x-\Delta x, t)}{(\Delta x)^2}

ここで、$ \Delta x$ は空間のステップサイズです。

時間微分の離散化

時間微分 $\frac{\partial u}{\partial t}$ を前進差分法で離散化します。これにより、次のようになります：

\frac{\partial u}{\partial t} \approx \frac{u(x, t+\Delta t) - u(x, t)}{\Delta t}

ここで、$\Delta t$ は時間のステップサイズです。

差分方程式の導出

これらの差分を熱伝導方程式に代入すると、次のような差分方程式が得られます：

\frac{u(x, t+\Delta t) - u(x, t)}{\Delta t} = \alpha \frac{u(x+\Delta x, t) - 2u(x, t) + u(x-\Delta x, t)}{(\Delta x)^2}

この式を $u(x, t+\Delta t)$ について解くと、次のようになります：

u(x, t+\Delta t) = u(x, t) + \alpha \frac{\Delta t}{(\Delta x)^2} \left[ u(x+\Delta x, t) - 2u(x, t) + u(x-\Delta x, t) \right]

実装された式

実際のコードでは、以下のように表されます：

new_temp[i] = temp[i] + kp * dt / dx**2 * (temp[i+1] - 2 * temp[i] + temp[i-1])

ここで、

temp[i] は時刻 $t$ における位置 $x$ での温度 $u(x, t)$ を表します。
new_temp[i] は時刻 $t+\Delta t$ における位置 $x$ での温度 $u(x, t+\Delta t)$ を表します。
kp は熱拡散率 $\alpha$ を表します。
dt は時間のステップサイズ $\Delta t$ を表します。
dx は空間のステップサイズ $\Delta x$ を表します。

この式により、次の時間ステップの温度分布を計算します。

数値解法と解析的解法を python で比較する例

数値解法では、有限差分法を用いて熱伝導方程式を離散化し、逐次的に温度分布を計算します。

google colab に実行例を置いています。

#!/usr/bin/env pythonimport numpy as np
import matplotlib.pyplot as plt
params = {'xtick.labelsize': 11, 'ytick.labelsize': 11, 'legend.fontsize': 8}
plt.rcParams['font.family'] = 'serif'

import imageio.v2 as imageio  # v2を明示的に使用
import os

def graph(x, y, x_label, y_label, title):
    """プロットを生成する関数"""
    plt.plot(x, y)
    plt.grid()
    plt.title(title)
    plt.xlabel(x_label)
    plt.ylabel(y_label)
    plt.show()

def initialize_temperature(imax, dx, xmax):
    """初期温度分布を生成する関数（方形波）"""
    temp = np.zeros(imax)
    for i in range(imax):
        x = dx * i
        temp[i] = 5 * (x / (xmax / 2)) if i < imax // 2 else 5 * (xmax - x) / (xmax / 2)
    return temp

def update_temperature(temp, imax, kp, dt, dx):
    """温度分布を更新する関数"""
    new_temp = temp.copy()
    for i in range(1, imax - 1):
        new_temp[i] = temp[i] + kp * dt / dx**2 * (temp[i+1] - 2 * temp[i] + temp[i-1])
    return new_temp

def analytical_solution(x, t, xmax, alpha, n_terms=1000):
    """解析的解を計算する関数"""
    L = xmax
    solution = np.zeros_like(x)
    for n in range(1, n_terms + 1):
        bn = 20 / (n * np.pi)**2 * (1 - (-1)**n)
        solution += bn * np.sin(n * np.pi * x / L) * np.exp(-alpha * (n * np.pi / L)**2 * t)
    return solution

def save_frame_with_analytical(X, temp, t, xmax, dt, kp, output_dir, frame_num):
    """各フレームを保存する関数"""
    analytical_temp = analytical_solution(X, t * dt, xmax, kp)
    plt.figure()
    plt.title(f"t = {t}")
    plt.plot(X, temp, "r", label="Numerical")
    plt.plot(X, analytical_temp, "b--", label="Analytical")
    plt.grid()
    plt.ylim(0,6)
    plt.xlabel("x")
    plt.ylabel("temperature (K)")
    plt.legend()
    plt.savefig(f"{output_dir}/frame_{frame_num:04d}.png")
    plt.show()
    plt.close()

def create_gif(output_dir, gif_name):
    """PNGファイルからGIFを生成する関数"""
    images = []
    for file_name in sorted(os.listdir(output_dir)):
        if file_name.endswith('.png'):
            file_path = os.path.join(output_dir, file_name)
            images.append(imageio.imread(file_path))
    imageio.mimsave(gif_name, images, duration=0.1)

# パラメータ設定
imax = 100
ntmax = 2500
dt = 0.004
kp = 1
xmax = 10
dx = xmax / imax
output_dir = 'frames'
gif_name = 'temperature_evolution.gif'

# 初期温度分布のプロット
X = np.linspace(0, xmax, imax)
temp = initialize_temperature(imax, dx, xmax)
graph(X, temp, "x", "temperature (K)", "Initial state")

# ディレクトリの作成
if not os.path.exists(output_dir):
    os.makedirs(output_dir)

# 温度分布のアニメーションフレームを保存
for t in range(ntmax):
    temp = update_temperature(temp, imax, kp, dt, dx)
    if t % 100 == 0:
        save_frame_with_analytical(X, temp, t, xmax, dt, kp, output_dir, t // 100)

# GIFを生成
create_gif(output_dir, gif_name)
print(f"GIF animation saved as {gif_name}")

アニメーションの生成

数値解と解析解を比較するために、各時刻における温度分布をアニメーションとして可視化します。上記のコードでは、各フレームをPNGファイルとして保存し、最後にそれらを結合してGIFアニメーションを生成しています。

初期条件がガウス分布の場合

初期条件がガウス分布の場合、フーリエ級数展開ではなく、ガウス関数の連続性を活用してフーリエ変換と逆フーリエ変換を用いて解析解を求めます。この方法は、ガウス関数がフーリエ変換を通じて再びガウス関数になるという特性を利用しています。以下に導出過程を示します。

初期条件のガウス分布

初期条件としてガウス分布を設定します：

u(x, 0) = \exp\left(-\frac{(x - \mu)^2}{2\sigma^2}\right)

ここで、$\mu$ はガウス分布の中心、$\sigma$ は標準偏差です。

ガウス分布のフーリエ変換

ガウス分布のフーリエ変換を考えます。フーリエ変換の定義に従うと、

\hat{u}(k, 0) = \int_{-\infty}^{\infty} u(x, 0) e^{-ikx} \, dx

ガウス分布 $u(x, 0)$ を代入します：

\hat{u}(k, 0) = \int_{-\infty}^{\infty} \exp\left(-\frac{(x - \mu)^2}{2\sigma^2}\right) e^{-ikx} \, d

ここで、変数変換 $ x' = x - \mu $ を行うと、

\hat{u}(k, 0) = \int_{-\infty}^{\infty} \exp\left(-\frac{x'^2}{2\sigma^2}\right) e^{-ik(x' + \mu)} \, dx'

これを整理すると、

\hat{u}(k, 0) = e^{-ik\mu} \int_{-\infty}^{\infty} \exp\left(-\frac{x'^2}{2\sigma^2}\right) e^{-ikx'} \, dx'

ガウス分布のフーリエ変換の特性により、

\int_{-\infty}^{\infty} \exp\left(-\frac{x'^2}{2\sigma^2}\right) e^{-ikx'} \, dx' = \sigma \sqrt{2\pi} \exp\left(-\frac{k^2 \sigma^2}{2}\right)

したがって、

\hat{u}(k, 0) = \sigma \sqrt{2\pi} \exp\left(-\frac{k^2 \sigma^2}{2}\right) e^{-ik\mu}

この結果からわかるように、ガウス分布のフーリエ変換には中心 $\mu$ が含まれています。ガウス分布の中心 $\mu$ はフーリエ変換の結果において複素指数関数の位相因子 $ e^{-ik\mu} $ として現れます。

波数空間での時間発展

次に、波数空間での時間発展を考えます。熱伝導方程式の波数空間での解は次のようになります：

\hat{u}(k, t) = \hat{u}(k, 0) \exp\left(-\alpha k^2 t\right)

これを先ほどのフーリエ変換した初期条件に適用すると、

\hat{u}(k, t) = \sigma \sqrt{2\pi} \exp\left(-\frac{k^2 \sigma^2}{2}\right) e^{-ik\mu} \exp\left(-\alpha k^2 t\right)

これを整理すると、

\hat{u}(k, t) = \sigma \sqrt{2\pi} \exp\left(-k^2 \left(\frac{\sigma^2}{2} + \alpha t\right)\right) e^{-ik\mu}

逆フーリエ変換

最後に、波数空間での解を逆フーリエ変換して元の空間に戻します。逆フーリエ変換を適用すると、得られる解は再びガウス分布の形になります：

u(x, t) = \frac{1}{\sqrt{2\pi (\sigma^2 + 2\alpha t)}} \exp\left(-\frac{(x - \mu)^2}{2 (\sigma^2 + 2\alpha t)}\right)

ここで、$\sigma$ は初期の標準偏差、$\alpha$ は熱拡散率、$\mu$ はガウス分布の中心です。

サンプルコード

下記がサンプルコードです。

import numpy as np
import matplotlib.pyplot as plt
import imageio.v2 as imageio  # v2を明示的に使用
import os

def graph(x, y, x_label, y_label, title):
    """プロットを生成する関数"""
    plt.plot(x, y)
    plt.grid()
    plt.title(title)
    plt.xlabel(x_label)
    plt.ylabel(y_label)
    plt.show()

def initialize_temperature_gaussian(imax, dx, xmax, mu, sigma):
    """初期温度分布を生成する関数（ガウス分布）"""
    temp = np.zeros(imax)
    for i in range(imax):
        x = dx * i
        temp[i] = np.exp(-0.5 * ((x - mu) / sigma) ** 2)
    return temp

def update_temperature(temp, imax, kp, dt, dx):
    """温度分布を更新する関数"""
    new_temp = temp.copy()
    for i in range(1, imax - 1):
        new_temp[i] = temp[i] + kp * dt / dx**2 * (temp[i+1] - 2 * temp[i] + temp[i-1])
    return new_temp

def analytical_solution_gaussian(x, t, alpha, mu, sigma):
    """ガウス分布の解析的解を計算する関数"""
    return np.exp(-0.5 * ((x - mu) / (sigma * np.sqrt(1 + 2 * alpha * t))) ** 2) / np.sqrt(1 + 2 * alpha * t)

def save_frame_with_analytical(X, temp, t, alpha, mu, sigma, output_dir, frame_num):
    """各フレームを保存する関数"""
    analytical_temp = analytical_solution_gaussian(X, t, alpha, mu, sigma)
    plt.figure()
    plt.title(f"t = {t}")
    plt.plot(X, temp, "r", label="Numerical")
    plt.plot(X, analytical_temp, "b--", label="Analytical")
    plt.grid()
    plt.ylim(0, 1.5)
    plt.xlabel("x")
    plt.ylabel("temperature")
    plt.legend()
    plt.savefig(f"{output_dir}/frame_{frame_num:04d}.png")
    plt.show()
    plt.close()

def create_gif(output_dir, gif_name):
    """PNGファイルからGIFを生成する関数"""
    images = []
    for file_name in sorted(os.listdir(output_dir)):
        if file_name.endswith('.png'):
            file_path = os.path.join(output_dir, file_name)
            images.append(imageio.imread(file_path))
    imageio.mimsave(gif_name, images, duration=0.1)

# パラメータ設定
imax = 100
ntmax = 2500
dt = 0.004
kp = 1
xmax = 10
dx = xmax / imax
mu = xmax / 2
sigma = 1
output_dir = 'frames_gaussian'
gif_name = 'temperature_evolution_gaussian.gif'

# 初期温度分布のプロット
X = np.linspace(0, xmax, imax)
temp = initialize_temperature_gaussian(imax, dx, xmax, mu, sigma)
graph(X, temp, "x", "temperature", "Initial state")

# ディレクトリの作成
if not os.path.exists(output_dir):
    os.makedirs(output_dir)

# 温度分布のアニメーションフレームを保存
for t in range(ntmax):
    temp = update_temperature(temp, imax, kp, dt, dx)
    if t % 100 == 0:
        save_frame_with_analytical(X, temp, t * dt, kp, mu, sigma, output_dir, t // 100)

# GIFを生成
create_gif(output_dir, gif_name)

print(f"GIF animation saved as {gif_name}")
import numpy as np
import matplotlib.pyplot as plt
import imageio.v2 as imageio  # v2を明示的に使用
import os

def graph(x, y, x_label, y_label, title):
    """プロットを生成する関数"""
    plt.plot(x, y)
    plt.grid()
    plt.title(title)
    plt.xlabel(x_label)
    plt.ylabel(y_label)
    plt.show()

def initialize_temperature_gaussian(imax, dx, xmax, mu, sigma):
    """初期温度分布を生成する関数（ガウス分布）"""
    temp = np.zeros(imax)
    for i in range(imax):
        x = dx * i
        temp[i] = np.exp(-0.5 * ((x - mu) / sigma) ** 2)
    return temp

def update_temperature(temp, imax, kp, dt, dx):
    """温度分布を更新する関数"""
    new_temp = temp.copy()
    for i in range(1, imax - 1):
        new_temp[i] = temp[i] + kp * dt / dx**2 * (temp[i+1] - 2 * temp[i] + temp[i-1])
    return new_temp

def analytical_solution_gaussian(x, t, alpha, mu, sigma):
    """ガウス分布の解析的解を計算する関数"""
    return np.exp(-0.5 * ((x - mu) / (sigma * np.sqrt(1 + 2 * alpha * t))) ** 2) / np.sqrt(1 + 2 * alpha * t)

def save_frame_with_analytical(X, temp, t, alpha, mu, sigma, output_dir, frame_num):
    """各フレームを保存する関数"""
    analytical_temp = analytical_solution_gaussian(X, t, alpha, mu, sigma)
    plt.figure()
    plt.title(f"t = {t}")
    plt.plot(X, temp, "r", label="Numerical")
    plt.plot(X, analytical_temp, "b--", label="Analytical")
    plt.grid()
    plt.ylim(0, 1.5)
    plt.xlabel("x")
    plt.ylabel("temperature")
    plt.legend()
    plt.savefig(f"{output_dir}/frame_{frame_num:04d}.png")
    plt.close()

def create_gif(output_dir, gif_name):
    """PNGファイルからGIFを生成する関数"""
    images = []
    for file_name in sorted(os.listdir(output_dir)):
        if file_name.endswith('.png'):
            file_path = os.path.join(output_dir, file_name)
            images.append(imageio.imread(file_path))
    imageio.mimsave(gif_name, images, duration=0.1)

# パラメータ設定
imax = 100
ntmax = 2500
dt = 0.004
kp = 1
xmax = 10
dx = xmax / imax
mu = xmax / 2
sigma = 1
output_dir = 'frames_gaussian'
gif_name = 'temperature_evolution_gaussian.gif'

# 初期温度分布のプロット
X = np.linspace(0, xmax, imax)
temp = initialize_temperature_gaussian(imax, dx, xmax, mu, sigma)
graph(X, temp, "x", "temperature", "Initial state")

# ディレクトリの作成
if not os.path.exists(output_dir):
    os.makedirs(output_dir)

# 温度分布のアニメーションフレームを保存
for t in range(ntmax):
    temp = update_temperature(temp, imax, kp, dt, dx)
    if t % 100 == 0:
        save_frame_with_analytical(X, temp, t * dt, kp, mu, sigma, output_dir, t // 100)

# GIFを生成
create_gif(output_dir, gif_name)

print(f"GIF animation saved as {gif_name}")

計算結果

まとめ

一次元の熱伝導方程式は数値計算とフーリエ級数展開やフーリエ変換のよい練習問題になりますね。

おまけ

熱伝導方程式はフランスのジョゼフ・フーリエ（Joseph Fourier）によって導入されました。フーリエは1822年に出版した著書「熱解析理論（Théorie analytique de la chaleur）」の中で、この方程式を提案し、熱の伝導現象を数学的に記述しました。フーリエの業績は、フーリエ級数やフーリエ変換、熱伝導の理論だけでなく、数学や物理学など様々な分野に広がっています。