9.3 Exercises
307
The initial condition is the famous and widely used Gaussian function with standard
deviation (or “width”) σ , which is here taken to be small, σ = 0.01, such that the
initial condition is a peak. This peak will then diffuse and become lower and wider.
Compute u(x, t) until u becomes approximately constant over the domain.
Filename: gaussian_diffusion.py.
Remarks Running the simulation with σ = 0.2 results in a constant solution
u ≈ 1 as t → ∞, while one might expect from “physics of diffusion” that the
solution should approach zero. The reason is that we apply Neumann conditions
as boundary conditions. One can then easily show that the area under the u curve
remains constant. Integrating the PDE gives
1
−1
∂u
∂t
dx = β
1
−1
∂ 2 u
∂x 2 dx .
Using the Gauss divergence theorem on the integral on the right-hand and moving
the time-derivative outside the integral on the left-hand side results in
∂
∂t
1
−1
u(x, t)dx = β
∂u
∂x
1
−1
= 0.
(Recall that ∂u/∂x = 0 at the end points.) The result means that
1
−1 udx remains
constant during the simulation. Giving the PDE an interpretation in terms of heat
conduction can easily explain the result: with Neumann conditions no heat can
escape from the domain so the initial heat will just be evenly distributed, but not leak
out, so the temperature cannot go to zero (or the scaled and translated temperature
u, to be precise). The area under the initial condition is 1, so with a sufficiently fine
mesh, u → 1, regardless of σ .
Exercise 9.7: Vectorize a Function for Computing the Area of a Polygon
Vectorize the implementation of the function for computing the area of a polygon
in Exercise 5.6. Make a test function that compares the scalar implementation in
Exercise 5.6 and the new vectorized implementation for the test cases used in
Exercise 5.6.
Hint Notice that the formula x 1 y 2 + x 2 y 3 + · · · + x n−1 y n =
n−1
i=0 x i y i+1 is
the dot product of two vectors, x[:-1] and y[1:], which can be computed
as numpy.dot(x[:-1], y[1:]), or more explicitly as numpy.sum(x[:-1]
*y[1:]).
Filename: polyarea_vec.py.
Exercise 9.8: Explore Symmetry
One can observe (and also mathematically prove) that the solution u(x, t) of the
problem in Exercise 9.6 is symmetric around x = 0: u(−x, t) = u(x, t). In such
a case, we can split the domain in two and compute u in only one half, [−1, 0]
or [0, 1]. At the symmetry line x = 0 we have the symmetry boundary condition
307
The initial condition is the famous and widely used Gaussian function with standard
deviation (or “width”) σ , which is here taken to be small, σ = 0.01, such that the
initial condition is a peak. This peak will then diffuse and become lower and wider.
Compute u(x, t) until u becomes approximately constant over the domain.
Filename: gaussian_diffusion.py.
Remarks Running the simulation with σ = 0.2 results in a constant solution
u ≈ 1 as t → ∞, while one might expect from “physics of diffusion” that the
solution should approach zero. The reason is that we apply Neumann conditions
as boundary conditions. One can then easily show that the area under the u curve
remains constant. Integrating the PDE gives
1
−1
∂u
∂t
dx = β
1
−1
∂ 2 u
∂x 2 dx .
Using the Gauss divergence theorem on the integral on the right-hand and moving
the time-derivative outside the integral on the left-hand side results in
∂
∂t
1
−1
u(x, t)dx = β
∂u
∂x
1
−1
= 0.
(Recall that ∂u/∂x = 0 at the end points.) The result means that
1
−1 udx remains
constant during the simulation. Giving the PDE an interpretation in terms of heat
conduction can easily explain the result: with Neumann conditions no heat can
escape from the domain so the initial heat will just be evenly distributed, but not leak
out, so the temperature cannot go to zero (or the scaled and translated temperature
u, to be precise). The area under the initial condition is 1, so with a sufficiently fine
mesh, u → 1, regardless of σ .
Exercise 9.7: Vectorize a Function for Computing the Area of a Polygon
Vectorize the implementation of the function for computing the area of a polygon
in Exercise 5.6. Make a test function that compares the scalar implementation in
Exercise 5.6 and the new vectorized implementation for the test cases used in
Exercise 5.6.
Hint Notice that the formula x 1 y 2 + x 2 y 3 + · · · + x n−1 y n =
n−1
i=0 x i y i+1 is
the dot product of two vectors, x[:-1] and y[1:], which can be computed
as numpy.dot(x[:-1], y[1:]), or more explicitly as numpy.sum(x[:-1]
*y[1:]).
Filename: polyarea_vec.py.
Exercise 9.8: Explore Symmetry
One can observe (and also mathematically prove) that the solution u(x, t) of the
problem in Exercise 9.6 is symmetric around x = 0: u(−x, t) = u(x, t). In such
a case, we can split the domain in two and compute u in only one half, [−1, 0]
or [0, 1]. At the symmetry line x = 0 we have the symmetry boundary condition