~poisson-eq-spherical-symmetry.py

import numpy as np
from matplotlib import pyplot as plt

dr = 0.01
f0 = 1.0
a = 10.0

r = np.arange(0, 50, dr)
u = np.where(r<=a, r**2/6 - a**2/2, -a**3/r/3) * f0
plt.xlabel('$r\'$', fontsize=12)
plt.ylabel('$u(r\')$', fontsize=12)
plt.plot(r, u)
plt.vlines(a, u.min(), u.max(), color='red', label='$r\'=a$')
plt.grid()
plt.legend()
plt.show()

~u-x-t-animation.py

import numpy as np
from matplotlib import pyplot as plt
from matplotlib import animation as animation

L = 1.0
c = 1.0
kappa = 0.2
dx = 0.001
dt = 0.01

x = np.arange(0, L+dx, dx)
us = []
ts = np.arange(0, 1.0+dt, dt)
for t in ts:
    # 波動方程式：固定端
    #u = np.sin(2*np.pi*x/L) * np.cos(2*np.pi*c*t/L)
    # 波動方程式：自由端
    #u = np.cos(2*np.pi*x/L) * np.cos(2*np.pi*c*t/L)
    # 拡散方程式：周囲の環境が一定
    #u = np.exp(-kappa * (2*np.pi/L)**2 * t) * np.sin(2*np.pi*x/L)
    # 拡散方程式：断熱
    u = np.exp(-kappa * (2*np.pi/L)**2 * t) * np.cos(2*np.pi*x/L)
    us.append(u)

def draw(us_, ts_, interval=50):
    u_min, u_max = np.min(us_), np.max(us_)
    def update(i):
        if i != 0:
            plt.cla()
        plt.title(r'$t = {:03f}$'.format(ts_[i]))
        plt.xlabel('$x$', fontsize=12)
        plt.ylabel('$u(x,t)$', fontsize=12)
        plt.ylim([u_min, u_max])
        plt.plot(x, us_[i], color='blue')
        plt.grid()
    fig = plt.figure()
    ani = animation.FuncAnimation(fig, update, interval = interval, frames = len(ts_))
    ani.save("wave-eq.gif", writer = 'imagemagick')
    plt.show()

draw(us, ts)

~wave-eq-introduction.py

import numpy as np
from matplotlib import pyplot as plt

x = np.arange(np.pi*13/32, np.pi*16/32, 0.001)
y = np.sin(x)
plt.plot(x, y, color='black')
i1, i2, i3 = len(x)//8, len(x)//2, len(x)*7//8
plt.plot(x[[i1, i2]], y[[i1, i2]], color='red')
plt.plot(x[[i1, i2, i2]], y[[i1, i1, i2]], color='red', lw=0.5)
plt.plot(x[[i2, i3]], y[[i2, i3]], color='blue')
plt.plot(x[[i2, i3, i3]], y[[i2, i2, i3]], color='blue', lw=0.5)
plt.text(x[i1]+0.01, y[i1]+0.0005, r'$\theta_{-}$', fontsize=12, color='red')
plt.text(x[i2]+0.03, y[i2]+0.0005, r'$\theta_{+}$', fontsize=12, color='blue')

plt.xticks([x[i1], x[i2], x[i3]], [r'$x-\Delta x$', r'$x$', r'$x+\Delta x$'])
plt.yticks([y[i1], y[i2], y[i3]], [r'$u(x-\Delta x,t)$', r'$u(x,t)$', r'$u(x+\Delta x,t)$'])

plt.grid()
plt.show()