波動方程式の有限差分法を半自動で導出(sympy)

sympyで差分化した2次元の波動方程式を定義して、式の整理を自動で実行する。

jupyterで実行時、最終的な実行結果は下記。

uの添え字は、1番目が時間、2番目と3番目がxとy。

左辺がu_{1,0,0}なので、次の時間ステップにおける波の振幅を表す。

右辺は、u{0,?,?}とu{-1,?,?}で構成されているため、既知の値から計算できる。

jupyter labから.mdで出力したものをそのままここに貼っても、texの数式がうまく表示されなかった。 下ページに従い、pythonコード形式に変換した。 ipynb ファイルを py ファイルに変換する方法 | ふじやまエッグのブリコルール日誌

#!/usr/bin/env python
# coding: utf-8

# # solve
# $$
# \frac{\partial^2 u}{\partial t^2} = c^2 \left[ \frac{\partial^2 u}{\partial x^2}+ \frac{\partial^2 u}{\partial y^2} \right]
# $$
# 
# with explicit finite difference method

# In[1]:


import sympy as sym


# ## define symbols

# discretized wave: $u_{time, x, y}$
# 
# mesh size of space: $ dx $
# 
# mesh size of time: $ dt $
# 
# speed of wave: $c$

# In[2]:


u =[ [[0 for i in range(3)] for j in range(3)]for k in range(3)]

for nt in range(3):
    for nx in range(3):
        for ny in range(3):
#             u[nt][nx][ny] = sym.Symbol("u[nt+("+str(nt-1)+"),nx+("+str(nx-1)+"),ny+("+str(ny-1)+")]") # for python code
            u[nt][nx][ny] = sym.Symbol("u_{"+str(nt-1)+","+str(nx-1)+","+str(ny-1)+"}") # for readable math expression
    pass
u[0][0][0]


# In[3]:


dx = sym.Symbol("dx")  # mesh size of space
dt = sym.Symbol("dt")  # mesh size of time
c = sym.Symbol("c")  # wave speed


# ## define partial derivative
# 
# $$ \frac{\partial^2 u}{\partial x^2} := \frac{-u_{t, x-1, y} + 2u_{t,x,y}-u_{t,x+1,y}}{dx^2} $$
# $$ \frac{\partial^2 u}{\partial y^2} := \frac{-u_{t, x, y-1} + 2u_{t,x,y}-u_{t,x,y+1}}{dx^2} $$
# $$ \frac{\partial^2 u}{\partial t^2} := \frac{-u_{t-1, x, y} + 2u_{t,x,y}-u_{t+1,x,y}}{dt^2} $$
# 

# In[4]:


d2x = sym.expand(((u[1][1][1] - u[1][0][1]) / dx - (u[1][2][1] - u[1][1][1]) / dx) / dx)
d2x


# In[5]:


d2y=sym.expand(((u[1][1][1] - u[1][1][0])/dx-(u[1][1][2] - u[1][1][1])/dx)/dx)
d2y


# In[6]:


d2t=sym.expand(((u[1][1][1] - u[0][1][1])/dt-(u[2][1][1] - u[1][1][1])/dt)/dt)
d2t


# ## define discretized wave equation

# In[7]:


eq=sym.Eq(d2t,c**2*(d2x+d2y))
eq


# ## solve the above eq. for $u_{1,0,0}$=u[2][1][1]

# In[8]:


rhs=sym.solve(eq,u[2][1][1])[0].expand().simplify()
rhs_simp=sym.collect(rhs,c).collect(dt).collect(dx)
rhs_simp


# # discretized wave eq. for computation
# 
# lhs is unknown value
# 
# rhs is constructed with known value

# In[9]:


sym.Eq(u[2][1][1],rhs_simp)