2次元の拡散方程式を陽解法で解く(Numpy+matplotlibによるアニメーション)

追記：数式が部分的に壊れるみたいです。見た目をどうにか整えるために hspace などちょっと汚いことをしてるんですが、Google Chart API がよくないみたい…

本記事の数式とコードの記述は主にRef.1(文末)に寄りかかっている。

2次元の拡散を表す偏微分方程式は次のように記述される。

$\frac{\partial u}{\partial t} = a \biggl( \frac{\partial^{2} u}{\partial x^{2}} + \frac{\partial^{2} u}{\partial y^{2}} \biggr)$

数値計算によってこの方程式の解を求めるためには離散化を行う。 $u(x, y, t)$ は $u_{i,j}^{(m)}$ と記述され、 $x=i\Delta x$ , $y=j\Delta y$ , $t = m \Delta t$ である(i, jによる微分を下付き文字によって簡単に記述している）。

uの離散的な一つの形は、時間に関する前進差分によって与えることができる。偏導関数をいきなり示すと、

f:id:xef:20130130145733p:plain

である。これを移項して時間前進的な表式を得る。

f:id:xef:20130130150156p:plain

安定性条件

誤差が時間成長しないことをフォン・ノイマンの安定性解析によって証明するそうですがよーわからん。とにかくタイムステップがこれより小さければ前進差分法は安定になるみたいです（詳しい話はリファレンスに投げっぱ）。

$\Delta t \leq \frac{1}{2a} \frac{(\Delta x \Delta y)^2}{(\Delta x)^2 + (\Delta y)^2}$

初期条件と境界条件

初期条件と境界条件を解きたいシチュエーションに合わせて随意に決めることができる。今回の初期条件は、 $r = (i \Delta x - 0.05)^{2} + (j \Delta y - 0.05)^{2}$ として、

$\begin{eqnarray} u_{i, j}^{(0)} = \left\{ \begin{array} 1: \hspace{10.0em} 0.05 \leq r \leq 0.1 \\ 0: \hspace{33.0em} otherwise \end{array} \right. \end{eqnarray}$

境界条件として、xとyの範囲を0～100（単位不明）の間に限って話を進める。すなわちxおよびyについて[0, 100]の範囲外で値は0とする。

Python+Numpy+matplotlibによる解曲面のプロット

# -*- coding: utf-8 -*-
import numpy as np
from matplotlib import pyplot as plt
from matplotlib import animation as animation
from pylab import *

# グローバル変数＞＜
dx = 0.01
dy = 0.01
a = 0.5

nx = int(1/dx)
ny = int(1/dy)

ui = np.zeros([nx, ny]) # initialize matrices
u = np.zeros([nx, ny])

def initialize_ui(condition, ui, dx, dy, nx, ny):
    for i in range(nx):
        for j in range(ny):
            if condition(i, j, dx, dy):
                ui[i, j] = 1

def evolve_ts(x, im, dx2, dy2, nx, ny, dt, a):
    global u, ui
    u[1:-1, 1:-1] = ui[1:-1, 1:-1] + a*dt*(
                        (ui[2:, 1:-1] - 2*ui[1:-1, 1:-1] + ui[:-2, 1:-1])/dx2
                      + (ui[1:-1, 2:] - 2*ui[1:-1, 1:-1] + ui[1:-1, :-2])/dy2
                    )
    ui = np.copy(u)
    im.set_array(ui)
    return im,

def main():
    global u, ui, dx, dy, nx, ny
    dx2 = dx ** 2
    dy2 = dy ** 2

    dt = dx2*dy2/(2*a*(dx2 + dy2))

    condition = (lambda i, j, dx, dy:
        (((i*dx-0.5)**2+(j*dy-0.5)**2<=0.1)
        & ( (i*dx-0.5)**2+(j*dy-0.5)**2>=.05))
    )


    initialize_ui(condition, ui, dx, dy, nx, ny)

    fig = plt.figure()
    img = plt.subplot(111)
    im = img.imshow(ui, cmap=cm.hot, interpolation='nearest', origin='lower')
    im.figure = fig
    fig.colorbar(im)
    ani = animation.FuncAnimation(fig, evolve_ts, 240,
            fargs=(im, dx2, dy2, nx, ny, dt, a,),
            blit=True, interval=50)
    # ani.save("diffusion.mp4", fps=15)
    plt.show()

if __name__ == "__main__":
    main()

Reference

主要参考文献というかトレース

Performance Python: Solving The 2D Diffusion Equation With numpy | t-square CC-BY-SA 3.0

ノイマンの安定性解析について

Von Neumann stability analysis - Wikipedia, the free encyclopedia

Methods of Numerical Simulation, Chapter 5 WikipediaのLax equivalence theoremのExternal Linksより

謙虚さが足りないんだよ。 2 次元の拡散方程式の数値解析シミュレーション（理論編）