In this article, we will make a math gif with Google Colaboratory. Let's solve the Coltveg-Dofries equation that describes a nonlinear wave by numerical calculation and see that a soliton appears.
With Google Colaboratory, you can work with Python in notebook format (as long as you have an internet connection) without installing anything on your PC.
In the past article [Let's make a math gif easily using Python-@wakabame] [1], there was a procedure for building an environment, but this is no longer necessary. The code is also executed on Google's cloud server, so you can use computational resources for free.
There are other features, so please refer to the official documentation for details. https://colab.research.google.com/notebooks/welcome.ipynb?hl=ja
This section is based on [Numerical calculation of KdV equation by pseudo-spectral method--without trying to do it] [2].
An initial value boundary value problem that describes the movement of shallow water wave height $ u
In numerical calculation, $ u_x $ and $ u_ {xxx} $ are obtained by Fourier transform using diff
of scipy.fftpack
in the spatial direction.
On the other hand, discretize in the time direction and use the solve_ivp
of scipy.integrate
to find the numerical solution at the next time step time by the implicit Runge-Kutta method.
import numpy as np
from numpy import pi, cos, linspace, arange
from numpy.fft import rfft, irfft
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from IPython.display import HTML
from scipy.integrate import solve_ivp
from scipy.fftpack import diff as psdiff
# Constant in Zabusky & Kruskal (1965)
DELTA = 0.022 ** 2
L = 2.0 # width of line domain
N = 256 # step width for spatial variable
T = 10 # last time
dt = 1.0e-2 # step width for time variable
# initial condition
x = linspace(0.0, L, N, endpoint=False)
u0 = cos(pi*x)
def kdv(t, u):
conv = u * psdiff(u, period=L)
disp = DELTA * psdiff(u, period=L, order=3)
dv = - conv - disp
return dv
t = np.arange(0, T, dt)
sol = solve_ivp(kdv, [0, T], u0, method='Radau', t_eval=t)
This section is based on [Visualizing Gradient Descent on Google Colaboratory-@yaju] [3].
The obtained solution sol.y
is a two-dimensional array arranged in the order of[space] [time]
, so transpose it so that it becomes ʻU [t] [x]`.
U = sol.y.T
frames = len(U)
When visualizing in a notebook, it is easy to see by checking from HTML
of ʻIPython.display`.
However, on Colaboratory, the connection is cut when a large amount of memory is consumed in one cell (?), So the gif frame is thinned out to fit in about 200.
%matplotlib inline
fig = plt.figure()
fig.set_dpi(100)
ax = fig.add_subplot(1,1,1)
def animate(t):
ax.clear()
plt.ylim([-1.2, 3.0])
p, = plt.plot(x, U[5*t])
anim = animation.FuncAnimation(fig, animate, frames=frames//5, interval=100, repeat=False)
HTML(anim.to_jshtml())
When exporting a gif file, the writer =" pillow "
option allows you to output it with almost no modification.
When the export is complete, open the tab on the left and you will find the gif file saved on the drive, so right click and download it
fig = plt.figure()
fig.set_dpi(100)
ax = fig.add_subplot(1,1,1)
def animate(t):
ax.clear()
plt.ylim([-1.2, 3.0])
p, = plt.plot(x, U[5*t])
anim = animation.FuncAnimation(fig, animate, frames=frames//5, interval=100, repeat=False)
anim.save("kdv_pseudo_spec_impl_rc.gif", writer="pillow")
This notebook is available at kdv.ipynb. If you have any comments, please do not hesitate to contact us. [1]:https://qiita.com/wakabame/items/c3648501eb0f2b921ddf [2]:https://iqujack-lequina.hatenablog.com/entry/2018/05/02/%E6%93%AC%E3%82%B9%E3%83%9A%E3%82%AF%E3%83%88%E3%83%AB%E6%B3%95%E3%81%AB%E3%82%88%E3%82%8BKdV%E6%96%B9%E7%A8%8B%E5%BC%8F%E3%81%AE%E6%95%B0%E5%80%A4%E8%A8%88%E7%AE%97 [3]:https://qiita.com/yaju/items/42d43d6d6cf6910e8701
Recommended Posts