7. Partial Differential Equations

7. Partial Differential Equations#

References:

Svein Linge & Hans Petter Langtangen: Programming for Computations – Python. Springer (2016).
- Chapter 5: Solving Partial Differential Equations

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

Partial derivative#

For a function with multiple inputs $f(x_1,...,x_n)$, a partial derivative is the derivative with respect to an input $x_i$ while other inputs are held constant. $$\frac{\partial f(x_1,...,x_n)}{\partial x_i} = \lim_{\Delta x \rightarrow 0} \frac{f(x_1,...,x_i+\Delta x,...,x_n) - f(x_1,...,x_i,...,x_n)}{\Delta x}$$

For example, for $$ f(x,y) = e^{x} \sin y, $$ partial derivatives are $$ \frac{\partial f(x,y)}{\partial x} = e^{x} \sin y $$ \frac{\partial f(x,y)}{\partial y} = e^{x} \cos y $$

x, y = np.meshgrid(np.linspace(0, 2, 20), np.linspace(-3, 3, 20))
f = np.exp(x) * np.sin(y)
dfdx = np.exp(x) * np.sin(y)
dfdy = np.exp(x) * np.cos(y)
# 3D plot
fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(projection='3d')
ax.plot_surface(x, y, f, cmap='viridis')
x0 = np.zeros_like(x)
x1 = np.ones_like(x)
ax.quiver(x, y, f, x1, x0, dfdx, color='g', length=0.2)
ax.quiver(x, y, f, x0, x1, dfdy, color='r', length=0.2)
plt.xlabel('x'); plt.ylabel('y');

_images/80ae0485adc2a4b0d34d76eb0b844dd5b2c5b0dfa3cba4235285a1ae83699d2e.png

Partial Differential Equation (PDE)#

A partial differential equation (PDE) is an equation that includes partial derivatives $\frac{\partial f(x_1,...,x_n)}{\partial x_i}$ of an unknown function $f(x_1,...,x_n)$.

A typical case is a function in space and time $y(x,t).$

A simple example is the diffusion equation (or heat equation) $$ \frac{\partial y(x,t)}{\partial t} = D \frac{\partial^2 y(x,t)}{\partial x^2} + g(y,x,t) $$ that describes the evolution of concentration (or temperature) $y$ in space $x$ and time $t$ with input $g(y,x,t)$ and the diffusion coefficient $D$.

Analytic Solution of PDE#

For a diffusion equation without input $$ \frac{\partial y}{\partial t} = D \frac{\partial^2 y}{\partial x^2}, $$ the solution is given by separation of variables.

By assuming that the solution is a product of temporal and spatial components $y(x,t) = u(t) v(x)$, the PDE becomes $$ \frac{\partial u(t)}{\partial t} v(x)= D u(t) \frac{\partial^2 v(x)}{\partial x^2} $$ \frac{1}{Du(t)}\frac{\partial u(t)}{\partial t} = \frac{1}{v(x)}\frac{\partial^2 v(x)}{\partial x^2} $$ For this equation to hold for any $t$ and $x$, a possible solution is for both sides to be a constant $C$. Then we have two separate ODEs: $$ \frac{du(t)}{dt} = C D u(t) $$ \frac{d^2v(x)}{dx^2} = C v(x) $$ for which we know analitic solutions.

By setting $C=-b^2\le 0$, we have $$ u(t) = C_0 e^{-b^2 Dt}, $$ v(x) = C_1 \sin bx + C_2 \cos bx. $$

Thus we have a solution $$ y(x,t) = e^{-b^2 Dt}( C_1 \sin bx + C_2 \cos bx) $$ where $b$, $C_1$ and $C_2$ are determined by the initial condition $y(x,0)$.
The equation tells us that higher spatial frequency components decay quiker.

Boundary condition#

For uniquely specifying a solution of a PDE in a bounded area, e.g., $x_0<x<x_1$, we need to specify either

the value $y(x,t)$ (Dirichlet boundary condition)
or the derivative $\frac{\partial y(x,t)}{\partial x}$ (Neumann boundary condition)

at the boundary $x_0$ and $x_1$ to uniquely determine the solution.

From PDE to ODE#

The standard way of dealing with space and time on a digital computer is to discretize them. For a PDE $$ \frac{\partial y}{\partial t} = D \frac{\partial^2 y}{\partial x^2} + g(y,x,t), $$ defined in a range $x_0 \le x \le x_0+L$, we consider a spatial discretization by $\Delta x = L/N$ x_i = x_0 + i\Delta x $(i=1,…,N)$ and $$ y_i(t) = y(x_i,t).$$

Then the PDE can be approximated by a set of ODEs $$ \frac{dy_i}{dt} = D \frac{y_{i+1}-2y_i+y_{i-1}}{\Delta x^2} + g(y_i,x_i,t). $$

Euler method#

First, let us solve the converted ODE by Eulre method.

def euler(f, y0, dt, n, args=()):
    """f: righthand side of ODE dy/dt=f(y,t)
        y0: initial condition y(0)=y0
        dt: time step
        n: iteratons
        args: parameter for f(y,t,*args)"""
    d = np.array([y0]).size  ## state dimension
    y = np.zeros((n+1, d))
    y[0] = y0
    t = 0
    for k in range(n):
        y[k+1] = y[k] + f(y[k], t, *args)*dt
        t = t + dt
    return y

def diff1D(y, t, x, D, inp=None):
    """1D Diffusion equaiton with constant boundary condition
    y: state vector
    t: time
    x: positions
    D: diffusion coefficient
    input: function(y,x,t)"""
    dx = x[1] - x[0]  # space step
    # shift to left and right and subtract
    d2ydx2 = (y[:-2] -2*y[1:-1] + y[2:])/dx**2
    # add 0 to both ends for Dirichlet boundary condition
    d2ydx2 = np.hstack((0, d2ydx2, 0))
    if inp == None:
        return D*d2ydx2
    else:
        return D*d2ydx2 + inp(y, x, t)

Lx = 5   # length
N = 20
x = np.linspace(0, Lx, N+1)
y0 = np.zeros_like(x)  # initial condition
y0[0] = 1   # 1 at left end

D = 0.1  # diffusion constant
dt = 0.1 # time step
ts = 15.   # time for solution
y = euler(diff1D, y0, dt, int(ts/dt), (x, D))  
plt.plot(x, y.T)
plt.xlabel("x"); plt.ylabel("y");

_images/e6e1ee3c8066570049c3bf6bf77f7b2f8fe20306ac6e0931ad16d468bb50ffda.png

plt.plot(y)
plt.xlabel("t"); plt.ylabel("y");

_images/6c561e3fd43e6bd5ad84509e90afeacd7475f1b1d64959a99dcf1ea47ab132d8.png

Let us try a finer spatial division

N = 50
dt = 0.1
x = np.linspace(0, Lx, N+1)
y0 = np.zeros_like(x)  # initial condition
y0[0] = 1   # 1 at left end
D = 0.1  # diffusion constant
y = euler(diff1D, y0, dt, int(ts/dt), (x, D))  
plt.plot(x, y.T)
plt.xlabel("x"); plt.ylabel("y");

_images/2c963c8357e238037817a5e9ca2308d56de22a3b4e47cdf18bf6e8efac51c86a.png

Stability#

For an accurate solution, $\Delta x$ have to be small. But making $\Delta x$ small can make the ODE stiff such that solution can be unstable.
It is known that stable solution by Euler method requires $$ \Delta t \le \frac{\Delta x^2}{2D}. $$

By odeint#

The above stability issue can be worked around by using an ODE library that implements automatic time step adaptation, such as odeint.

from scipy.integrate import odeint

Lx = 5   # length
N = 50
x = np.linspace(0, Lx, N+1)
y0 = np.zeros_like(x)  # initial condition
y0[0] = 1
D = 0.1  # diffusion constant
dt = 0.1
t = np.arange(0, ts, dt)
y = odeint(diff1D, y0, t, (x, D))
p = plt.plot(x, y.T)
plt.xlabel("x"); plt.ylabel("y");

_images/97a2480fefecf8e3e7c6bfc7881718ed3698006d0060ee9ee51a18415ec02c43.png

Let us see the case with an input.

def pulse(y, x, t):
    """1 for 1<x<2 at 1<t<2"""
    return (1<x)*(x<2)*(1<t)*(t<2)*1.

N = 50
x = np.linspace(0, Lx, N+1)
y0 = np.zeros_like(x)  # initial condition
ts = 10
t = np.arange(0, ts, dt)
y = odeint(diff1D, y0, t, (x, 1, pulse))
p = plt.plot(x, y.T)
plt.xlabel("x"); plt.ylabel("y");

_images/4d1d3906fb0a7e95d3f2284530d964da8dfc256e6e288e202e5173831876f342.png

# plot of each point in time
p = plt.plot(t, y)
plt.xlabel("t"); plt.ylabel("y");

_images/05c8fdb0db3f18856437070d26c51efe34488fc3d06122911f3bda99db0408f0.png

# plot in space and time
plt.imshow(y.T, origin='lower', extent=(0, ts, 0, Lx))
plt.xlabel("t"); plt.ylabel("x");

_images/81afa0d2d956bba9da285210a761b28237b89b6d7a0e55bda10b797f764974f8.png

3D plot and animation#

%matplotlib notebook

# plot in 3D
fig = plt.figure()
ax = fig.add_subplot(projection='3d')
T, X = np.meshgrid(t, x)
ax.plot_surface(T, X, y[:,:N+1].T, cmap='viridis')
plt.xlabel("t"); plt.ylabel("x"); 

from matplotlib import animation

fig = plt.figure()
frames = []  # prepare frame
for i, ti in enumerate(t):
    p = plt.plot(x, y[i])
    plt.xlabel("x"); plt.ylabel("y");
    frames.append(p)    
anim = animation.ArtistAnimation(fig, frames, interval = 10)

Boundary Conditions#

As briefly mentioned above, there are two typical ways of specifying the boundary condition:

Dirichlet boundary condition: the value $y(x_0,t)=y_0$ or $y(x_N,t)=y_N$
- concentration at an open end of a tube to a water bath.
- voltage of an electric cable attached to the ground or a strong battery.
Neumann boundary condition: the derivative $\left.\frac{\partial y(x,t)}{\partial x}\right|_{x_0}=y'_0$ or $\left.\frac{\partial y(x,t)}{\partial x}\right|_{x_N}=y'_N$
- concentration at a closed end of a tube (no molecular flow).
- voltage at a dead end of an electric cable.

We already implemented Dirichlet boundary condition in diff1D().

Let us implement a Neumann boundary condition at the right end $\left.\frac{\partial y(x,t)}{\partial x}\right|_{x_N}=0$ and see the difference.

A simple way is to use a one-sided differentiation $$ \frac{\partial y}{\partial x} \simeq \frac{y_{N+1}-y_{N}}{\Delta x}=0$$ which reduces to $y_{N+1}=y_{N}$ and $$ \frac{dy_N}{dt} = D \frac{y_{N-1} - y_N}{\Delta x^2} + g(y,x,t).$$

A better way is to use a symmetric differentiation $$ \frac{\partial y}{\partial x} \simeq \frac{y_{N+1}-y_{N-1}}{2\Delta x}=0$$ which reduces to $y_{N+1}=y_{N-1}$ and $$ \frac{dy_N}{dt} = 2D \frac{y_{N-1} - y_N}{\Delta x^2} + g(y,x,t).$$

def diff1DN(y, t, x, D, inp=None):
    """1D Diffusion equaiton with boundary conditions
    y_0=0 on the left and dy/dx|x_N=0 on the right end.
    y: state vector
    t: time
    x: positions
    D: diffusion coefficient
    input: function(y,x,t)"""
    # spatial step
    dx = x[1] - x[0]
    # shift array to left and right
    d2ydx2 = (y[:-2] -2*y[1:-1] + y[2:])/dx**2
    # Dirichlet on the left, Neumann on the right end
    d2ydx2 = np.hstack((0, d2ydx2, 2*(y[-2] - y[-1])/dx**2))
    if inp == None:
        return(D*d2ydx2)
    else:
        return(D*d2ydx2 + inp(y, x, t))

N = 100
x = np.linspace(0, 1, N+1)
y0 = np.zeros_like(x)  # initial condition
y0 = np.sin(np.pi*x)
t = np.arange(0, 1, 0.01)
y = odeint(diff1DN, y0, t, (x, 1))
plt.figure()
plt.plot(x, y.T)
plt.xlabel("x"); plt.ylabel("y");

Wave Equation#

PDE with a second-order dynamics can represent traveling waves. The wave equation has the standard form $$ \frac{\partial^2 u}{\partial t^2} + d\frac{\partial u}{\partial t} = c^2 \frac{\partial^2 u}{\partial x^2}$$ where $c$ is the wave speed and $d$ is the decay rate.

We convert the second-order system to a set of first order systems by consiering a vector $y=(u,v)$ representing the amplitude $u$ and its change rate $v$: $$ \frac{\partial u}{\partial t} = v $$ \frac{\partial v}{\partial t} = c^2 \frac{\partial^2 u}{\partial x^2} - dv $$

def wave1D(y, t, x, c, d, inp=None):
    """1D wave equaiton with constant boundary
    y: state vector hstack(u, v)
    t: time
    x: positions
    c: wave speed
    input: function(y,x,t)"""
    n = int(len(y)/2)
    u, v = y[:n], y[n:]
    dx = x[1] - x[0]
    # finite different approximation
    d2udx2 = (u[:-2] -2*u[1:-1] + u[2:])/dx**2
    d2udx2 = np.hstack((0, d2udx2, 0))  # add 0 to both ends
    if inp == None:
        return np.hstack((v, c**2*d2udx2 - d*v))
    else:
        return np.hstack((v, c**2*d2udx2 - d*v + inp(y, x, t)))

Lx = 10
N = 100
x = np.linspace(0, Lx, N+1)
y0 = np.zeros(2*(N+1))  # initial condition
c = 1.  # wave speed
d = 0.1
ts = 30
dt = 0.1
t = np.arange(0, ts, dt)
y = odeint(wave1D, y0, t, (x, c, d, pulse))
plt.figure()
plt.plot(x, y[:,:N+1].T)  # only first dimension
plt.xlabel("x"); plt.ylabel("u"); 

# plot in 3D
fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(projection='3d')
T, X = np.meshgrid(t, x)
ax.plot_surface(T, X, y[:,:N+1].T, cmap='viridis')
plt.xlabel("t"); plt.ylabel("x"); 

fig = plt.figure()
frames = []  # prepare frame
for i, ti in enumerate(t):
    p = plt.plot(x, y[i,:N+1])
    plt.xlabel("x"); plt.ylabel("y");
    frames.append(p)    
anim = animation.ArtistAnimation(fig, frames, interval = 10)

/opt/anaconda3/lib/python3.12/site-packages/matplotlib/animation.py:892: UserWarning: Animation was deleted without rendering anything. This is most likely not intended. To prevent deletion, assign the Animation to a variable, e.g. `anim`, that exists until you output the Animation using `plt.show()` or `anim.save()`.
  warnings.warn(

Reaction-Diffusion Equation#

Linear waves can decay, reflect, and overlay.
Some nonlinear waves, like neural spike, can travel without decaying. FitzHugh-Nagumo model can be embedded in a diffusive axon model $$ \frac{\partial v}{\partial t} = v - \frac{v^3}{3} - w + D\frac{\partial^2 v}{\partial x^2} + I $$ \frac{\partial w}{\partial t} = \phi (v + a - bw) $$

def fhnaxon(y, t, x, inp, D=1, a=0.7, b=0.8, phi=0.08):
    """FitzHugh-Nagumo axon model
    y: state vector hstack(v, w)
    t: time
    x: positions
    I: input current
    D: diffusion coefficient
    """
    n = int(len(y)/2)
    v, w = y[:n], y[n:]
    # finite difference approximation
    d2vdx2 = (v[:-2] -2*v[1:-1] + v[2:])/(x[1] - x[0])**2
    d2vdx2 = np.hstack((0, d2vdx2, 0))  # add 0 to both ends
    dvdt = v - v**3/3 - w + D*d2vdx2 + inp(y, x, t)
    dwdt = phi*(v + a -b*w)
    return(np.hstack((dvdt, dwdt)))

def pulse2(y, x, t):
    """1 for 5<x<10 at 5<t<10"""
    return (5<x)*(x<10)*(5<t)*(t<10)

Lx = 50
N = 100
x = np.linspace(0, Lx, N+1)
y0 = np.zeros(2*(N+1))  # initial condition
y0[0:N+1] = -2  # reset 
ts = 50
dt = 0.5
t = np.arange(0, ts, dt)
y = odeint(fhnaxon, y0, t, (x, pulse2))
plt.figure()
p = plt.plot(x, y[:,:N+1].T) # plot in space
plt.xlabel("x"); plt.ylabel("v");

plt.figure()
plt.imshow(y[:,:N+1].T, origin='lower', extent=(0, ts, 0, Lx))
plt.xlabel("t"); plt.ylabel("x");

# plot in 3D
fig = plt.figure()
ax = fig.add_subplot(projection='3d')
T, X = np.meshgrid(t, x)
ax.plot_surface(T, X, y[:,:N+1].T, cmap='viridis')
plt.xlabel("t"); plt.ylabel("x"); 

fig = plt.figure()
frames = []  # prepare frame
for i, ti in enumerate(t):
    p = plt.plot(x, y[i,:N+1])
    plt.xlabel("x"); plt.ylabel("y");
    frames.append(p)    
anim = animation.ArtistAnimation(fig, frames, interval = 10)

PDE in 2D space#

We can also model dynamics in two or higher dimension spate.

The diffusion equation of a variable $y(x_1,x_2,t)$ in the 2D space is given by: $$\frac{\partial y}{\partial t} = D\frac{\partial^2 y}{\partial x_1^2} + D\frac{\partial^2 y}{\partial x_2^2} + g(y, x, t)$$

Using the Laplace operator (or Laplacian): $$\nabla^2 = \frac{\partial^2}{\partial x_1^2} + \frac{\partial^2}{\partial x_2^2}$$ the diffusion equation is represented as $$\frac{\partial y}{\partial t} = D\nabla^2 y + g(y, x, t)$$

Diffusion equation in 2D space#

def diff2D(y, t, x, D, inp=0):
    """1D Diffusion equaiton with constant boundary condition
    y: states at (x0, x1)
    t: time
    x: positions (x0, x1): transposed meshgrid
    D: diffusion coefficient
    inp: function(y,x,t) or number"""
    n = x.shape[:2]  # grid size (n0, n1)
    nn = n[0]*n[1]  # grid points n0*n1
    y = y.reshape(n)  # linear to 2D
    dx = [x[1,0,0]-x[0,0,0], x[0,1,1]-x[0,0,1]]  # space step
    # Laplacian with cyclic boundary
    Ly = (np.roll(y,-1,0) + np.roll(y,1,0) + np.roll(y,-1,1) + np.roll(y,1,1) - 4*y)/dx[0]**2
    Ly[0,:] = Ly[-1:0] = Ly[:,0] = Ly[:,-1] = 0  # Dirichlet boundary
    dydt = D*Ly + (inp(y,x,t) if callable(inp) else inp)
    return dydt.ravel()

def pulse3(y, x, t):
    """1 for 0.2<x0<0.3 at 1<t<2"""
    return (0.2<x[0])*(x[0]<0.3)*(0.4<x[1])*(x[1]<0.6)*(1<t)*(t<2)

L = 10
N = 20
X = np.meshgrid(np.linspace(0,L,N), np.linspace(0,L,N))
x = np.array(X).transpose(2,1,0)  # (x0,x1) in last dimention
y0 = np.zeros((N,N))  # initial condition
y0[2:4,5:8] = 1  # 
ts = 5
dt = 0.1
t = np.arange(0, ts, dt)

y = odeint(diff2D, y0.ravel(), t, (x, 1))
plt.figure(figsize=(8,3))
for i in range(10):
    plt.subplot(2, 5, i+1)
    plt.imshow(y[5*i].reshape(N,N).T)
    plt.title('t={:1.1f}'.format(t[5*i]))
plt.tight_layout()

fig = plt.figure()
frames = []  # prepare frame
for i, ti in enumerate(t):
    p = plt.imshow(y[i].reshape(N,N).T)
    frames.append([p])    
anim = animation.ArtistAnimation(fig, frames, interval = 10)
plt.xlabel("x1"); plt.ylabel("x2");

fig = plt.figure(figsize=(6,6))
ax = fig.add_subplot(projection='3d')
frames = []  # prepare frame
for i, ti in enumerate(t):
    p = ax.plot_surface(X[0], X[1], y[i].reshape(N,N).T, cmap='viridis')
    frames.append([p])    
anim = animation.ArtistAnimation(fig, frames, interval = 10)
plt.xlabel("x1"); plt.ylabel("x2");

PDE libraries#

Solution of PDE is computationally intensive, especially in 2D or higher dimensional space. For practical computation, it is better to use specialized libraries for solving PDEs, such as: