熱方程式の数値解析(3) - らんだむな記憶をpythonで再実装。
わりと計算が速かった。octaveは遅いのかな...。
#! /usr/bin/env python # -*- coding: utf-8 -*- import numpy as np import matplotlib.pyplot as plt def initial_func(xx): length_xx = len(xx) y = np.zeros(length_xx) for i in range(length_xx): x = xx[i] if x <= -1 or x >= 1: y[i] = 0 else: y[i] = np.exp(- 1 / (1 - x**2)) return y def laplacian(uu): length_uu = len(uu) y = np.zeros(length_uu) for i in range(length_uu): u = uu[i] left_u = 0 right_u = 0 if i - 1 >= 0: left_u = uu[i - 1] if i + 1 < length_uu: right_u = uu[i + 1] y[i] = (left_u + right_u - 2*u) / 2; return y class HeatEquation(object): def __init__(self): self.w = 10 # half width of bar self.dx = 0.05 self.dt = self.dx**2 self.T = 10 # total elapsed time [sec] self.n = int(self.T/self.dt) # step num # sampling on x-axis self.xx = np.linspace(-self.w, self.w, int(2*self.w/self.dx)) self.u = None self.__setup() def solve(self): length_xx = len(self.xx) # solve heat equation by forward Euler method for step in range(self.n-1): t = self.dt * step # current time self.u[step + 1, :] = self.u[step, :] + laplacian(self.u[step, :]) def show(self): plt.plot(self.xx, self.u[0,:], "r") plt.plot(self.xx, self.u[int(self.n/5),:], "b") plt.plot(self.xx, self.u[int(2*self.n/5),:], "k") plt.plot(self.xx, self.u[int(3*self.n/5),:], "g") plt.plot(self.xx, self.u[int(4*self.n/5),:], "y") plt.plot(self.xx, self.u[self.n-1,:], "m") plt.axis([-10, 10, 0, 0.4]) plt.show() def __setup(self): length_xx = len(self.xx) # solutions at each time self.u = np.zeros([self.n, length_xx]) # initial value self.u[0,:] = initial_func(self.xx) ################################################################################ heat = HeatEquation() heat.solve() heat.show()
結果の描画:
―――――・・・
この辺の実践的な説明と理論的な説明はともに金子先生の本にある。
実践的説明
理論的説明
波動方程式の数値解析 - らんだむな記憶にぽそっと書いた条件$\frac{(\Delta x)^2}{2 \Delta t} = D^\prime \le \frac{1}{2}$は「Courant-Friedrichs-Lewyの条件」と呼ばれるもので、1928年には理論的に与えられたものらしい。
http://gdz.sub.uni-goettingen.de/dms/load/img/?PPN=GDZPPN002272636&IDDOC=29349←たぶんこれ。ドイツ語なので厳しい...。
Lewyってどっかで聞いたことあるなと思ったけど、Lewy's example - Wikipedia, the free encyclopediaで聞く名前だな。一般に微分方程式は初期値の微小の摂動に対して解の変化も微小であることを期待したい。そういうのが期待できないと、排水管とかで水量がちょっと変動すると排水管に思わぬ圧力がかかって破裂したりするかもしれないから。解の存在と一意性に加え、この安定性をもってHadamardのwell-posednessという概念になる、はず。
共著者のRichard CourantはCourant-HilbertのCourantのようだ。Kurt Otto FriedrichsはMollifier(軟化子)やFriedrichs extensionの人だろう。自身がないので一応調べた。後者はFriedrichs extension in nLabに書いてある。Hilbert空間の閉対称作用素の自己共役実現を保証する大事な定理だ。とりあえず、大物3人の共同論文ということが分かったが、とんでもないな...。
このCourant-Friedrichs-Lewyの条件は金子先生の本によると、差分スキームに依存するもののようで、上記のような前進差分スキームだと制限が出てくるが、後退差分スキーム($u(t,x) - u(t - \Delta t,x)$)だと$D^\prime$の値は何であっても解は安定するらしい。なんてこった...。必ずしも方程式に固有のものではないということのようだ。ランダムウォークから熱方程式を導く場合には$D^\prime = \frac{1}{2}$として結構自然な形で現れるので、何か絶対的な条件なのかなという気もしたのだが。
同書には、浅水波の伝播を記述するKorteweg–de Vries方程式の数値計算についても書いてあって、孤立波解(ソリトン)の様子を見ることができるようだ。差分スキームの都合上、初期値の座標だけでなく、速度に相当するような情報も与えないとならないようだ。気が向いたら実装してみたい気もする。
mKdV方程式との架け橋であるMiura変換はRobert M. Miura氏が導いたのだな。やっぱし日系(日系3世)なのか。ふむ。
