第7回
１次元波動方程式の数値解法
（2021/03/03 最終修正）
Embedded Files
ここでは，１次元波動方程式(wave equation)をコンピュータを用いて解くために差分化を行い，コンピュータを用いて波動方程式を解くアルゴリズムを示す．前回の拡散方程式に比べると、時間微分が2階になっているため、多少ややこしいが考え方は同じである。特に、初期条件の与え方がわかりにくいかも知れない。
さて，天下り的だが1次元波動方程式は弦の振動を（おおよそ）あらわす．次の動画は、1次元波動方程式をここで紹介する数値解法で解いた結果をアニメーションとしたものだ。 （こちらは、C言語＋あるグラフィックライブラリで作成したアニメーション）
下記方程式 (W) における u(x,t) は弦の（静止位置からの）変位をあらわす．ここでは1次元波動方程式を数値計算を用いて解き，その解 u(x,t) をアニメーションとして表示する事を目標とする．つまり，コンピュータ内部に弦の振動を再現し，それをアニメーションとして表示するわけだ．
さて，次の1次元波動方程式の初期値境界値問題、
但し、c >0, f(0)=0, f(L)=0, g(0)=0, g(L)=0．
を考える。境界条件，初期条件が適切に与えられれば (W) の解が弦の振動を表現していると言えるだろう．さて，数値計算する為にまず(W)を離散化する必要がある。これはこれまでの常微分方程式の場合、また拡散方程式と共通で，要は微分を差分で置き換える必要がある．分割方法は以下の２種類が考えられるが，
ここでは時間方向の離散化は (1) を，空間方向の離散化は (2) を採用します．（前回の拡散方程式と同じ）
[0,l] 区間を N 等分するとし，上記離散化方法に従い，h = l/N, τ>0 とし，xj = (j - 1/2)h, tk = kτ として，
と書くことにすると、(W)の第一式は次のように離散化される。
λ=cτ/h として、まとめると
という式が得られる。2階微分の差分化については，このページ最後のスライドを参照．この式は、時刻 tk-1, tkでの U から次の時刻 tk+1での U が計算できるという形をしている。したがって、初期条件として、 Uj0,Uj1 が必要となる。まず Uj0は(W)の第三式より、
となる。 Uj1 を決める方法はいくつか考えられるが、例えば次のようにする。まず、テイラーの公式を用いると、
となる。(W)より
であり，さらに t=0 でも(W)の第一式が成り立つとすると、
となる。この式の右辺は２階差分商を用いて
と近似できる。つまり，
と求まる．
境界条件は（W）の第二式だが，まずは境界条件について触れる必要がある．ここでは，境界の処理の為に，U0k,UN+1k という架空のセルを用意することにする．
様々な境界条件があるが，今回の（W) における境界条件は Dirichlet 条件（この場合各時間 0 という値に固定されている．一般の場合，例えば a という値に固定されている場合には，次の式の右辺 0 が a となる．）と呼ばれるもの．空間離散化の方法として（１）を採用した事を思い出すと，
となればよいことが分かる．すなわち，
となる。（この辺り、拡散方程式と同じ）
まとめると、以下のようになる．
但し、 λ≦1 が安定性のための必要十分条件（安定性条件とは，数値計算がうまく行く為の条件）.
めでたく差分方程式が得られた。(2)をコンピューターを用いて解くのがここでの目標である。状況を、次の絵を使って説明する。
縦軸方向を時間の進行方向とし，横軸方向を空間方向とする．赤色の円盤で示した位置の値 Ukj が，他の4点，Uk-1j ,Uk-1j-1 ,Uk-1j+1 ,Uk-2j ,から決定される様子．領域の端以外では上記のイメージで良いが，端においては境界条件の考慮が必要となる（前述）．
1次元波動方程式を解く為のアルゴリズムは次のようになる。なお、安定性条件を満たすため、λをパラメータとして先に与えて、τは与えられたλから求めるようにする。
アルゴリズム中の uj, old_uj, new_uj はそれぞれ実数型の配列として定義すると良い．境界条件処理用に配列要素はN個よりも余分に２つ必要となり，計N+2個必要となる事に注意（new_uj については，実際には境界処理用の要素は用いない）．例えば uj について，Python風に配列を定義すると，
u = [0.0]*(N+2)
となる．このようにすれば，u[0], u[1], ... , u[N+1] が用意される．実際そのように用意すれば良い．
以下のアルゴリズムのイメージ．new_u という配列要素を u および old_u という配列要素から求める．
初期条件における、new_u, u, old_u の関係は、old_u が時刻0の弦の形状、すなわち f(x)であり与えられているもの、u が時刻 τ における弦の形状であり、f(x)およびg(x)から決定されるもの、new_u はそれらから求まる時刻 2τ における弦の形状であり未知であるものとなっている。
アルゴリズム
アルゴリズム
（１）初期パラメータ設定
N, T, λ, c, l を設定する (λ ≦ 1)
h = l/N, τ=(λh)/c, M = T/τ
x, new_u, u, old_u を配列（リスト）として用意
（２） 初期値設定
j := 1,....,N の順に
x[j] = (j - 0.5)h
old_u[j] = f(x[j])
を繰り返す
old_u[0]=-old_u[1], old_u[N+1]=-old_u[N]（境界条件：初期値決定のため）
j := 1,....,N の順に
u[j]=old_u[j]+τg(x[j])+(λ*λ)/2*(old_u[j-1]-2old_u[j]+old_u[j+1])
を繰り返す
（３）
k := 0,1,.....,M-1 の順に
(old_u[1] ～ old_u[N] を画面に表示)
u[0]=-u[1], u[N+1]=-u[N]（境界条件）
j := 1,....,N の順に
new_u[j]=2u[j]-old_u[j]+(λ*λ)*(u[j-1]-2u[j]+u[j+1])
を繰り返す
j := 1,....,N の順に
old_u[j] = u[j]
を繰り返す
j := 1,....,N の順に
u[j] = new_u[j]
を繰り返す
を繰り返す
課題1（演習）
課題1（演習）
l=1, c=1, f(x)=2*x*(1-x), g(x)=0, u(0,t)=u(1,t)=0 として (2) を解くプログラムを作成せよ。(Nは100程度)
Code7-1.py : 1次元波動方程式のシミュレーション
Code7-1.py : 1次元波動方程式のシミュレーション
# Code7-1.py
import math
import matplotlib.pyplot as plt
import matplotlib.animation as animation
T = 5.0
M = 5000
tau = T/float(M)
L = 1.0
N = 100
h = L/float(N)
C = 1.0
Lambda = C*tau/h
INTV = 100
print(Lambda)
def f(x):
return 2*x*(1.0 - x)
#return math.exp(-100*(x-0.4)*(x-0.4))
def g(x):
return 0.0
#return -200*(x - 0.4)*f(x)
x, old_u, u, new_u = [0.0]*(N+2), [0.0]*(N+2), [0.0]*(N+2), [0.0]*(N+2)
fig = plt.figure()
ims = []
for j in range(1, N+1):
x[j] = (j - 0.5)*h
old_u[j] = f(x[j])
old_u[0] = -old_u[1]
old_u[N+1] = -old_u[N]
for j in range(1, N+1):
u[j] = old_u[j] + tau*g(x[j]) + (Lambda*Lambda/2.0)*(old_u[j-1] - 2*old_u[j] + old_u[j+1])
for k in range(0, M+1):
t = k*tau
if k%INTV == 0:
im = plt.plot(x[1:N+1], u[1:N+1], linewidth=2, color="blue")
ims.append(im)
u[0] = -u[1]
u[N+1] = -u[N]
for j in range(1, N+1):
new_u[j] = 2*u[j] - old_u[j] + (Lambda*Lambda)*(u[j-1] - 2*u[j] + u[j+1])
for j in range(1, N+1):
old_u[j] = u[j]
for j in range(1, N+1):
u[j] = new_u[j]
plt.title("7-1")
plt.xlabel("x")
plt.ylabel("u(x,t)")
ani = animation.ArtistAnimation(fig, ims, interval=100, repeat=False)
plt.show()
#ani.save("anim.gif", writer="pillow", fps=10)
課題2
課題2
l=1, c=1, u(0,t)=u(1,t)=0 はそのままで、f(x) および g(x) を他のもにかえてみよ。（もっと局在化した初期値からはじめると、波の伝播が観察される。）
f(x)=exp(-100*(x-0.4)*(x-0.4))
とし、異なる g(x) を初期条件としたものを、
(a)waveeq2
(b)waveeq3
この時の g(x) は(a),(b)それぞれどのようなものか考えよ。
実際の弦の振動現象と照らし合わせ，f(x) および g(x) の意味を再考せよ．
課題3
課題3
弦の振動ではなく，膜の振動を考えたい．すなわち，１次元波動方程式ではなく２次元波動方程式を考える．どのような方程式になるであろうか？
2次元波動方程式を数値計算で解き、その結果を OpenGL という可視化ライブラリを用いて可視化したサンプルは以下の様。
これに近いものを、おまけ回に載せたので参考にして欲しい。
差分近似について．
差分近似について．
Page updated
Google Sites
Report abuse