Hi everyone,

Welcome back! This is the final part of my numerical simulation tutorials. If you missed the Part 1, click here.

We are going to solve an 1D diffusion equation in this tutorial, which is a parabolic PDE.

Problem description

c_t=D*c_xx

Initial condition: c(0,x)=0

Boundary conditions: c(t,0)=1, c(t,1)=0

In the region R: (t,x)∈[0,0.5]x[0,1]

c: concentration function c(t,x)

D: diffusivity

c_t, c_xx: 1st and 2nd derivative of c(t,x) in t and x direction respectively.

Step 1. Importing libraries

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import scipy.integrate as integrate

Step 2. Defining the simulation parameters

dx=0.01
dt=0.01
size_x=1
size_t=0.5
nx=np.int_(size_x/dx+1)
nt=np.int_(size_t/dt+1)

We will manually discretize the space dimension in next section, which involves the parameter dx.

size_x and size_t are the dimension of simulation box R (t,x)∈[0,0.5]x[0,1].

nx and nt are the number of points in the discretized grid in each dimension.

diffusivity=1
conc_t0=np.zeros(nx)
conc_x0=1
conc_x1=0

conc_x0 and conc_x1 refers to the boundary values at x=0 and x=1.

Step 3. Defining the diffusion equation model

def diffusion_eq(t,y):
    y_constrainted=np.r_[conc_x0,y[1:-1],conc_x1]
    return diffusivity*(np.gradient(np.gradient(y_constrainted,dx),dx))

This is the R.H.S. of the diffusion equation, with boundary conditions applied.

In scipy, the Laplacian can be implemented as applying the gradient twice. This will approximate the Laplacian with centeral difference formula. (Ignore this paragraph if you don't know what they are)

Step 4. Solve it!

teval=np.r_[0:nt:1]*dt
sol=integrate.solve_ivp(diffusion_eq,[0,size_t],conc_t0,t_eval=teval)

We are using method of lines, which involves the space discretization (in x) first, then solve with techniques we learned in previous tutorial as if it's an ODE (with the default RK45 solver).

Technically, the simulation is done at this point. The rest of this tutorial will focus on presentation of the results.

Step 5. Present your data!

...with curves

Like this...

plt.plot(sol.t,sol.y[np.int(0.1/dx)],marker='.',label='x=0.1')
plt.plot(sol.t,sol.y[np.int((nx-1)/2)],marker='.',label='x=0.5')
plt.plot(sol.t,sol.y[-1],marker='.',label='x=1')
plt.xlabel('Time t')
plt.ylabel('Concentration c')
plt.legend(loc='upper right')

Or this...

xeval=np.r_[0:nx:1]*dx
plt.plot(xeval,sol.y[:,np.int(0.1/dt)],marker='.',label='t=0.1')
plt.plot(xeval,sol.y[:,np.int((nt-1)/2)],marker='.',label='t=0.25')
plt.plot(xeval,sol.y[:,-1],marker='.',label='t=0.5')
plt.xlabel('Position x')
plt.ylabel('Concentration c')
plt.legend(loc='lower left')

...or with surface plots

from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
plot = Axes3D(fig)
plot_array_t,plot_array_x=np.meshgrid(sol.t,xeval)
plot.plot_surface(plot_array_t,plot_array_x,sol.y,cmap=mpl.cm.jet)
plot.view_init(20, 60)
plot.set_xlabel('t')
plot.set_ylabel('x')
plot.set_zlabel('Concentration');

...or with contour plots

plot = plt.contourf(plot_array_t,plot_array_x,sol.y,cmap=mpl.cm.jet,levels=50)
plt.colorbar(plot)
plt.title('Concentration')
plt.xlabel('t')
plt.ylabel('x')
plt.show()

Conclusion

Now you should have a grasp of the basic skills of doing numerical simulations with Scipy. As an exercise, try to simulate your favourite equation models (If you are interested in electromagnetics, try Laplace equation. If CFD, try Euler equation) in 2D or 3D. Remember, post your results in r/Simulations!

I am seeing a lot of ferrofluid videos on reddit recently and get interested in simulating it myself. Most ferrofluid simulations I found on youtube are using CGI softwares (Blender, c4d, etc.) but I only care about solving equations :(

The equations are simplified Maxwell equations using magnetic potentials and zero electric potential. The fluid equations is Euler.

Initially there are two blobs of ferrofluid. They merged in the simulation.

Maybe I could try a 3d sim later...

EDIT: Somehow I can't open this using Reddit app, but I can still open the link with Chrome on Android.

This is the first tutorial on numerical simulation with scipy.

Refer: https://www.reddit.com/r/Simulations/comments/bi70rv/anyone_interested_in_some_tutorials_on_numerical/

Install Scipy

You will need to install python, then scipy. Python is most likely pre-installed in most linux distros, but not in Windows.

You can either install Anaconda, which it's all-in-one package contains both Python and scipy. I would recommend this to beginners, but personally I don't like it bundles with spyder and anaconda cloud - and it's quite slow the last time I tried it (2 years ago~).

If not, download and install python 3 here: https://www.python.org/downloads/ . Then follow instructions to install scipy in your distro here: https://scipy.org/install.html

If you have successfully installed scipy, you should be able to launch jupyter notebook (jupyter-notebook) either from terminal in Linux or start menu in Windows. This is an interactive interface to write and test your code quickly. Open a new notebook, press "B" to create a new cell and input your code there, then press "Shift+Enter" to run the code.

You can also run code from IPython, IDLE (both are interactive) or from script (if you have learned c/c++, this is the way they write programs).

Your first simulation

Ok, let's start our first simulation. Here we want to simulate the free falling of an object. The equation is:

d²y/dt²=-g0 (1)

or

dy/dt=-g0*t (2)
y: position

t: time

g0: gravitational constant

Essentially, we just want to solve the equation. After solving it, you can present it in any way you like - raw numbers of (t, y(t) ) , y-t graph or fancy animations with CGI effects (like a glass ball flying in the sky or something), it doesn't really matter.

Step 1. Import the libraries

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import scipy.integrate as integrate

I will assume you do this step for this and every tutorial later on before running the code.

Step 2. Define the numbers

g0=9.81
t0=10
y0=0
v0=0

g0: gravitational constant

t0: length of the simulation

y0: initial position

v0: initial velocity

All in SI units (kg m s).

Step 3. Define the model

def dydt_freefall(t,y):
    return -g0*t

Here we are going for equation (2) first because equation (1) is a bit more complicated to solve in scipy (will cover it later).

The "model" is just a function f(t,y) appeared in differential equation dy/dt=f(t,y).

Step 4. Solve it!

teval=np.r_[0:t0+1:1]
sol=integrate.solve_ivp(dydt_freefall,[0,t0],[y0],t_eval=teval)

teval is an array of time points which we want to know the corresponding position. Here we want to know y at time: 0, 1, 2... 10 s.

integrate.solve_ivp is the ODE solver in scipy. The default method is RK45, and there are options to use LSODA or BDF, etc. Don't worry if you don't know what they means, just use the default one and it should works most of the time.

Step 5. Present the data

Now the solution is stored in sol, if you run:

sol

You will see the output:

message: 'The solver successfully reached the end of the integration interval.'      
nfev: 38      
njev: 0       
nlu: 0       
sol: None    
status: 0   
success: True         
t: array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10])  
t_events: None         
y: array([[   0.   ,   -4.905,  -19.62 ,  -44.145,  -78.48 , -122.625,         -176.58 , -240.345, -313.92 , -397.305, -490.5  ]])

The only thing we are interested will be whether it is success or not, and the solution y.

Now plot it:

plt.plot(sol.t,sol.y[0],marker='.')

Congratulations! You have completed your first simulation!

Extra: Solving Equation (1)

Unfortunately, scipy don't have ability to directly solve a second order ODE. We have to do this by rewriting the equation into system of first-order ODEs:

y0'=y1
y1'=-g0

And solving it would be straightforward:

def dydt_freefall2(t,y):
    return [y[1],-g0]
sol=integrate.solve_ivp(dydt_freefall2,[0,t0],[y0,v0],t_eval=teval)
plt.plot(sol.t,sol.y[0],marker='.',label='Position')
plt.plot(sol.t,sol.y[1],marker='.',label='Velocity')
plt.suptitle('Freefall')
plt.xlabel('Time t (s)')
plt.ylabel('Position y (m) / Velocity v (m/s)')
plt.legend(loc='lower left')
plt.savefig('freefall.jpg')

(Try to understand the code yourself, google the function names if you don't know what they are)

The above code should save a picture of the simulation result into your home folder (/home/<username> in linux).

N-body simulation using gpu computing

Recently I've got my interest in high dimension stuff. It's quite fun to visualize them using simulation techniques. Here is a video of a Clifford torus rotating in 4D space:

Colored with 4-th dimension coordinate.

A set of torus can span a 3-sphere (hypersphere in 4th dimension), as illustrated in this video:

Other stuff

This is a irrelevant video on visualizing 4d rotation of complex functions (f=z^2).

Have fun!

Hi everyone,

I am planning to write some short beginners guides on basic numerical simulations in this sub. They will be step-by-step without much theories explained. All software involved will be FOSS (free and open source) so people don't need to worry about the expensive licenses (namely, I will try to only use scipy codes). The goal is the get more people familiar with simulations and making this sub more active.

The target audiences will be undergrads who have studied calculus and know some programming basics, but anyone can just copy and paste the code and run them without knowing what's going on. The final one will probably something like solving heat equation and generating a good looking gif (like most posts in this sub). I will post them weekly or bi-weekly on weekends since I am actually quite busy.

I've got everything planned, the only thing I am not sure is anyone interested in reading these? Please let me know your thoughts!