Diffusion equation animation with NumPy

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

Diffusion.gif

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

Recommended Posts

Diffusion equation animation with NumPy
Animation with matplotlib
Animation with matplotlib
Moving average with numpy
Getting Started with Numpy
Learn with Cheminformatics NumPy
Matrix concatenation with Numpy
Hamming code with numpy
Regression analysis with NumPy
Extend NumPy with Rust
Kernel regression with Numpy only
I wrote GP with numpy
CNN implementation with just numpy
Artificial data generation with numpy
[Python] Calculation method with numpy
Equation of motion with sympy
Debt repayment simulation with numpy
Implemented SMO with Python + NumPy
Stick strings together with Numpy
Handle numpy arrays with f2py
Use OpenBLAS with numpy, scipy
Bubble sort with fluffy animation
Python3 | Getting Started with numpy
Implementing logistic regression with NumPy
Create plot animation with Python + Matplotlib
Perform least squares fitting with numpy.
Draw a beautiful circle with numpy
Implement Keras LSTM feedforward with numpy
Extract multiple elements with Numpy array
Easy animation with matplotlib (mp4, gif)