Numerical methods for a nonlinear parabolic problem
Numerical methods for a nonlinear parabolic problem
Jeffrey M. Connors
February 20, 2015
The problem is to find u(x, t) : R2 → R solving
ut = uxx + f(u), (t > 0), (1)
f(u) = u(u − β)(1 − u), 0 < β <
1
2
, (2)
u(−∞, t) = 1, (3)
u(∞, t) = 0, (4)
u(x, 0) = u0(x), (traveling wave). (5)
We shall describe two numerical methods to estimate the solution. We must truncate the domain for
computational purposes, so we estimate the boundary conditions as u(−L, t) = 1 and u(L, t) = 0.
Here, L 1 is an adjustable parameter (that should be set near the top of the code). Let M be the
number of interior nodes xj in space, j = 1, 2, . . . ,M, with uniform grid spacing Δx = 2L/(M +1);
x0 = −L and xM+1 = L. Here, the goal is for u(xj, t) ≈ uj(t). We discretize (1) in space as follows:
∂
∂t
uj(t) =
1
Δx2 (uj+1 − 2uj + uj−1) + f(uj), (6)
for j = 1, …M; u0 = 1 and uM+1 = 0. We expect that when Δx is small, |u(xj, t)−uj(t)| ≈ CΔx2
for some C > 0.
Consider first using a time-stepping method such as Crank-Nicolson (below). Let Δt > 0 be
the size of the time step, with N total time steps to be taken. Let our approximations now be
unj
≈ u(xj, tn) with u0j
= u0(xj) and define tn+1/2 = tn +Δt/2. The method is derived by noting
that
ut(tn+1/2) =
u(tn+1) − u(tn)
Δt
+ O
Δt2
and then estimating all other terms in (6) at time tn+1/2. Let un+1/2
j = (un+1
j + unj
)/2. There are
a number of ways to define the method, but it is convenient to use
un+1
j − unj
Δt
=
1
Δx2
un+1/2
j+1 − 2un+1/2
j + un+1/2
j−1
+ f(un+1/2
j ), (7)
for j = 1, 2, . . . ,M. Note that these are only the interior nodes; boundary values get absorbed into
the right-hand side of the equations as data.
We are not done yet. We need to address the nonlinearity. Two approaches will be taken for
this purpose, as outlined below.
Extrapolated Crank-Nicolson
We replace the nonlinear term f(un+1/2
j ) in (7) with the second-order accurate extrapolation
E(f)n+1/2
j =
3
2
f
unj
−
1
2
f
un−1
j
.
1
The result is
un+1
j − unj
Δt
=
1
Δx2
un+1/2
j+1 − 2un+1/2
j + un+1/2
j−1
+ E(f)n+1/2
j ,
which we manipulate into the form
un+1
j −
Δt
2Δx2
un+1
j+1 − 2un+1
j + un+1
j−1
= unj
+
Δt
2Δx2
unj
+1 − 2unj
+ unj
−1
+ΔtE(f)n+1/2
j . (8)
This is implemented by solving a linear system at each time step of the form
(I + γD)un+1 = Fn = (I − γD)un +ΔtE (f)n+1/2 + 2γ[1 0 . . . 0]T
with I the M ×M identity matrix, D a diffusion matrix with γ = Δt/(2Δx2), uk a column vector
with components ukj
and E (f)n+1/2 a column vector with components E(f)n+1/2
j . The last term
on the right-hand side is due to boundary conditions.
Newton’s Method
In this case, we iterate to resolve the nonlinearity using Newton’s method to find Φ(un+1) = 0,
where
Φ(z) = (I + γD)z − (I − γD)un − ΔtG((z + un)/2) − 2γ[1 0 . . . 0]T . (9)
Here z has components zj , j = 1, …M and G(η) = [f(η1) . . . f(ηM)]T . The iterates will be z(k),
k = 0, 1, . . .. We will choose the initial guess z(0) = un at each time step. Then apply the usual
Newton iteration:
DΦ(z(k))Δz(k) = −Φ(z(k))
z(k+1) = z(k) +Δz(k).
One can use (2) and some standard calculus to show that
DΦ(z(k))Δz(k) =
I + γD − ΔtG(z(k))
Δz(k),
G(z(k)) = diag (g1, g2, . . . , gM) ,
gj =
1
8
−3
z(k)
j + unj
2
+ 4(1 + β)
z(k)
j + unj
− 4β
, j = 1, . . . ,M.
Each iteration requires one to solve a linear system. It is preferable to eliminate Δz(k) and solve
directly for z(k+1). This is accomplished via the following equation:
I + γD − ΔtG(z(k))
z(k+1) =(I − γD)un +ΔtG((z(k) + un)/2)
− ΔtG(z(k))z(k) + 2γ[1 0 . . . 0]T .
(10)
One iterates until z(k+1) − z(k) < , where > 0 is a user-defined parameter in the code. In
rare instances, it may be necessary to check (in addition) that Φ(z(k+1)) < . Finally, Newton’s
method may not converge in general, so a maximum k-value, say KMAX, should be specified such
that the code exits upon reaching k > KMAX.
2
