1

Poisson’s equation

Suppose that we want to describe the distribution of heat throughout a region Ω. Let h(x) represent
the temperature on the boundary of Ω (∂Ω), and let g(x) represent the initial heat distribution at
time t = 0. If we let f (x, t) represent any heat sources/sinks in Ω, then the flow of heat can be
described by the boundary value problem (BVP)

ut = 4u + f (x, t), x ∈ Ω, t > 0,

u(x, t) = h(x), x ∈ ∂Ω, (1.1)
u(x, 0) = g(x).

When the source term f does not depend on time, there is often a steady-state heat distribution u∞
that is approached as t → ∞. This steady state u∞ is a solution of the BVP

4u + f (x) = 0, x ∈ Ω,
(1.2)
u(x, t) = h(x), x ∈ ∂Ω.

This last partial differential equation, 4u = −f , is called Poisson’s equation. This equation
is satisfied by the steady-state solutions of many other evolutionary processes. Poisson’s equation is
often used in electrostatics, image processing, surface reconstruction, computational fluid dynamics,
and other areas.

Poisson’s equation in two dimensions

Consider Poisson’s equation together with Dirichlet boundary conditions on a square domain R =
[a, b] × [c, d]:

uxx + uyy = f, x in R ⊂ R2 ,
(1.3)
u = g, x on ∂R.

Let a = x0 , x1 , . . . , xn = b be a partition of [a, b], and let c = y0 , y1 , . . . , yn = d be a partition of [c, d].

Suppose that there are n + 1 evenly spaced points, so that n is the number of subintervals in each
dimension, and xi , yj are given by

xi = a + ih,
yj = c + jh,

1
2 Lab 1. Poisson’s equation

for i, j = 0, . . . , n − 1, where h = xi − xi−1 = yi − yi−1 . We look for an approximation Ui, j on the

grid {(xi , yj )}ni,j=0 .
Recall that

4u = uxx (x, y) + uyy (x, y)

u(x + h, y) − 2u(x, y) + u(x − h, y)
=
h2
u(x, y + h) − 2u(x, y) + u(x, y − h)
+ + O(h2 ).
h2
We replace 4 with the finite difference operator 4h , defined by
Ui+1, j − 2Ui, j + Ui−1, j Ui, j+1 − 2Ui, j + Ui, j−1
4h Uij = 2
+ ,
h h2
1
= 2 (Ui−1, j + Ui+1, j + Ui, j−1 + Ui, j+1 − 4Ui, j ).
h
Then the set of equations

4h Uij = fij , i, j = 1, . . . , n − 1,

can be written in matrix form as

AU + p + q = f.
A is a block tridiagonal matrix, given by
 
T I
I T I 
1  .. .. ..

(1.4)
 
h2 
 . . . 

 I T I
I T

where I is the (n − 1) × (n − 1) identity matrix, and T is the (n − 1) × (n − 1) tridiagonal matrix

−4 1
 
 1 −4 1 
 

 . .. . .. . .. .

 
 1 −4 1 
1 −4

The vector U is given by

 1   
U U1, j
 U2 
 where U j =  U2, j
 
U =  for each j, 1 ≤ j ≤ n − 1.
   
U n−1 Un−1, j

The vectors p and q come from the boundary conditions of (1.3), and are given by
 1   1 
p q
 ...   ... 
p= 
, q = 
 
,

pn−1 q n−1
 3

where  
g0, j
 0 
1  
pj = 2  ...  , 1 ≤ j ≤ n − 1,
 
h  
 0 
gn, j
and      
g1,0 g1,n 0
 g2,0  g2,n
 0
1 
..
 1  
..
 
0
q = 2 , q n−1
= 2 , q j =  ...  , 2 ≤ j ≤ n − 2.
     
h

 .  h

 .
  
gn−2,0  gn−2,n  0
gn−1,0 gn−1,n 0

Problem 1. Complete the function poisson_square by implementing the finite difference

method
AU + p + q = f
where the matrix A and the vectors p and q are defined following (1.4). Use your function
to solve the BVP
∆u = 0, x ∈ [0, 1] × [0, 1],
3
(1.5)
u(x, y) = x , (x, y) ∈ ∂([0, 1] × [0, 1]).

Plot the 3D solution with n = 100.

Poisson’s equation and conservative forces

In physics Poisson’s equation is used to describe the scalar potential of a conservative force. In
general
∆V = −f
where V is the scalar potential of the force, or the potential energy a particle would have at that
point, and f is a source term. Examples of conservative forces include Newton’s Law of Gravity
(where matter become the source term) and Coulomb’s Law, which gives the force between two
charge particles (where charge is the source term).
In electrostatics the electric potential is also known as the voltage, and is denoted by V. From
Maxwell’s equations it can be shown that that the voltage obeys Poisson’s equation with the electric
charge density (like a continuous cloud of electrons) being the source term:
ρ
∆V = − ,
ε0
where ρ is the charge density and ε0 is the permissivity of free space, which is a constant that we’ll
leave as 1.
Usually a non zero V at a point will cause a charged particle to move to a lower potential,
changing ρ and the solution to V . However, in this analysis we’ll assume that the charges are fixed
in place.
4 Lab 1. Poisson’s equation

Figure 1.1: The solution of (1.5).

Suppose we have 3 nested pipes. The outer pipe is attached to "ground," which usually we
define to be V = 0, and the inner two have opposite relative charges. Physically the two inner pipes
would function like a capacitor.
The following code will plot the charge distribution of this setup.

import matplotlib.colors as mcolors

def source(X,Y):
"""
Takes arbitrary arrays of coordinates X and Y and returns an array of the ←-
same shape
representing the charge density of nested charged squares
"""
src = np.zeros(X.shape)
src[ np.logical_or(
np.logical_and( np.logical_or(abs(X-1.5) < .1,abs(X+1.5) < .1) ,abs(Y) ←-
< 1.6),
np.logical_and( np.logical_or(abs(Y-1.5) < .1,abs(Y+1.5) < .1) ,abs(X) ←-
< 1.6))] = 1
src[ np.logical_or(
 5

np.logical_and( np.logical_or(abs(X-0.9) < .1,abs(X+0.9) < .1) ,abs(Y) ←-

< 1.0),
np.logical_and( np.logical_or(abs(Y-0.9) < .1,abs(Y+0.9) < .1) ,abs(X) ←-
< 1.0))] = -1
return src

#Generate a color dictionary for use with LinearSegmentedColormap

#that places red and blue at the min and max values of data
#and white when data is zero

def genDict(data):
zero = 1/(1 - np.max(data)/np.min(data))
cdict = {'red': [(0.0, 1.0, 1.0),
(zero, 1.0, 1.0),
(1.0, 0.0, 0.0)],
'green': [(0.0, 0.0, 0.0),
(zero, 1.0, 1.0),
(1.0, 0.0, 0.0)],
'blue': [(0.0, 0.0, 0.0),
(zero, 1.0, 1.0),
(1.0, 1.0, 1.0)]}
return cdict

a1 = -2.
b1 = 2.
c1 = -2.
d1 = 2.
n =100
X = np.linspace(a1,b1,n)
Y = np.linspace(c1,d1,n)
X,Y = np.meshgrid(X,Y)

plt.imshow(source(X,Y),cmap = mcolors.LinearSegmentedColormap('cmap', genDict(←-

source(X,Y))))
plt.colorbar(label="Relative Charge")
plt.show()

The function genDict scales the color values to be white when the charge density is zero. This
is mostly to help visualize where there are neutrally charged zones by forcing them to be white. You
may find it useful to also apply it when you solve for the electric potential.
With this definition of the charge density, we can solve Poisson’s equation for the potential
field.
6 Lab 1. Poisson’s equation

0 1.0
0.8
20 0.6
0.4
0.2

Relative Charge
40
0.0
60 0.2
0.4
80 0.6
0.8
1.0
0 20 40 60 80

Figure 1.2: The charge density of the 3 nested pipes.

Problem 2. Solve
∆V = −ρ(x, y), x ∈ [−2, 2] × [−2, 2],
(1.6)
V (x, y) = 0, (x, y) ∈ ∂([−2, 2] × [−2, 2]).

for the electric potential V. Use the source function (−ρ) defined above. Plot the 2D solution
using n = 100.
 7

Figure 1.3: The electric potential of the 3 nested pipes.