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.

1. Use Google Colaboratory

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

2. Numerical calculation of the Coltweg-Dofries equation

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 , called the Coltveg-Dofries equation $ \begin{cases} &u_t + \alpha u u_x + \beta u_{xxx} = 0, \quad & t>0, x \in (0, L),\\ &u(t,0) = u(t,L),\quad &t>0,\\ &u(0,x) = u_0(x), \quad &x \in (0,L) \end{cases} $$ Think about. Here, $ L $, $ \ alpha $, $ \ beta $ are positive constants, and the parameters when Zabsky Kruskal discovers a solitary wave $ L = 2.0 $, $ \ alpha = 1.0 $, $ \ beta = Assume that it is 0.022 ^ 2 $ and the initial value is given by $ u_0 (x) = \ cos (\ pi x) $.

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 required libraries

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

Parameter and initial value settings

# 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)

Calculation of $ u_t $ using Fourier transform

def kdv(t, u):
    conv = u * psdiff(u, period=L)
    disp = DELTA * psdiff(u, period=L, order=3)
    dv = - conv - disp
    return dv

Calculation of numerical solutions by the implicit Runge-Kutta method

t = np.arange(0, T, dt)
sol = solve_ivp(kdv, [0, T], u0, method='Radau', t_eval=t)

3. Let's visualize

This section is based on [Visualizing Gradient Descent on Google Colaboratory-@yaju] [3].

Preparation of numerical solution

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)

Visualization with notebook

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())

Export gif file

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

[PYTHON] Let's easily make a math gif using Google Colaboratory