things to do

I wanted to simulate a differential equation including a Laplacian with python, but the other day I came up with a method to process the difference of the Laplacian other than the for statement (reinvented the wheel), so I wrote an animation of the diffusion equation. It was.

When I looked it up again after I came up with it, there was a person who was upward compatible with what I did a few years ago. The article was very helpful: Speeding up numerical calculations using NumPy / SciPy: Application 1

The differential equation to be solved is a one-dimensional diffusion equation expressed by $ \frac{\partial }{\partial t}\psi(x,t)= D\frac{\partial^2}{\partial x^2}\psi(x,t) $ The initial condition was given by Gaussian, and the time evolution was Euler's method for the time being.

The code I wrote

`Diffusion_1D_ani.ipynb`


%matplotlib nbagg
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.animation import ArtistAnimation

# calculate Laplacian with periodic boundary condition
def Lap(L, M):
    # expand vector
    VL = [L[0]]
    VR = [L[-1]]
    Le = np.hstack([VR ,L, VL])
    Lc = np.dot(M, Le) # central diff
    return Lc[1:N+1]

# generate 1-D lattice and initialize
N = 64
x = np.linspace(-10, 10, N)
psi = np.exp(-x**2/7)

# time and diffusion const
dt = 0.05
T = np.arange(0.0, 37, dt)
D = 3

# make (N+2)x(N+2) matrix to calc Laplacian
I = np.eye(N+2,N+2)
S = np.vstack((I[1:],np.array([0]*(N+2))))
M = S + S.T -2*I 

# make animation
fig = plt.figure()
ims = []
for i in range(len(T)):
    psi += dt*(D*Lap(psi, M))
    line = plt.plot(x, psi, 'b')
    ims.append(line)

ani = ArtistAnimation(fig, ims, interval=5, blit=True, repeat=True)

plt.show()

Execution example

Ingenuity

When calculating the Laplacian, the values at the right and left ends of the original $ N $ dimensional array $ \ psi $ were attached to the left and right ends with np.hstack to form a $ N + 2 $ dimensional array. As a result, the diffusion term including the periodic boundary condition could be calculated automatically by taking the inner product with the $ (N + 2) \ times (N + 2) $ matrix $ M $ prepared in advance. The part added for the boundary condition is unnecessary, so it is discarded at the end of the function.

Summary

Problem: I wanted to make it visible for the time being, so I did not consider the numerical stability and error of 1 mm. I think there is something in NumPy that can do numerical integration quickly and safely, so I'll investigate and do something about it. I will add to these issues when any progress is made in the future.

Future: I think I'd like to simulate a reaction-diffusion system because it seems that a 2D version can be created by extending this 1D version appropriately. (If it is not an equation that includes a very ugly nonlinear term, it will be somehow ...)

Postscript:

(11/23) If we compromised on the forward difference, we could put a non-linear term in 3 lines. It seems that a one-dimensional Navier-Stokes equation can be created.

def nl(L, N):
    Lf = np.hstack((L,L[0]))
    return L * np.diff(Lf)

In the case of forward difference, the calculation diverged immediately if not only the Reynolds number but also the initial conditions and the step size of time were not carefully selected. What should i do

Diffusion equation animation with NumPy