Solving Partial Differential Equations

We all know Theano as a forefront library for Deep Learning research. However, it should be noted that Theano is a general purpose numerical computing library, like Numpy. Hence, in this post, we will look at the implementation of PDE simulation in Theano.

The Laplace Equation

We will look at a simple PDE example, the Laplace Equation:

$\nabla^2 \phi = 0$

In other words, this is an second order PDE, as, recall that Laplacian in calculus is the divergence of gradient of a function:

$\nabla \cdot \nabla \phi = 0$

particularly, in two dimension:

$\frac{\partial^2 \phi}{\partial x^2} + \frac{\partial^2 \phi}{\partial y^2} = 0$

This simple equation could be solved by using Finite Difference scheme [1].

Note that in this example, we are ignoring the boundary value problem.

Solving Laplace Equation in Numpy

The Finite Difference solution of Laplace Equation is to repeatedly averaging the neighbors of a particular point:

$\phi_{i, j} = \frac{1}{4} \left( \phi_{i+1, j} + \phi_{i+1, j} + \phi_{i, j+1} + \phi_{i, j-1} \right)$

This iterative solution is very simple to implement in Numpy. But first, let’s give this problem an initial condition:

$\phi = -5 \, \exp (x^2 + y^2)$

Visualizing the function:

Initial

Now, let’s implement this. First, we create the mesh, the solution space:

# Create 21x21 mesh grid
m = 21
mesh_range = np.arange(-1, 1, 2/(m-1))
x, y = np.meshgrid(mesh_range, mesh_range)

# Initial condition
U = np.exp(-5 * (x**2 + y**2))

We then create an indexing scheme that select the point north, west, south, and east of any given point. We do this so that we could implement this in vectorized manner.

References

Mitra, Ambar K. “Finite difference method for the solution of Laplace equation.” preprint.