AMATH 586 Hai Zhu Student ID: 1323888

Final Project

Hai Zhu

due date: June 9, 2014

In past homeworks, you have studied the ODE

1
v 0 (t) = g(v) (Ex0.0a)

where
g(v) = v(α − v)(v − 1) (Ex0.0b)
with 0 < α < 1.
We can turn (Ex0.0a) into a PDE in one space dimension and time by letting v(x, t) vary in
space and adding spatial diffusion, obtaining
1
vt (x, t) = κvxx (x, t) + g(v). (Ex0.0c)

This is a scalar reaction-diffusion equation, a variant of the Cahn-Allen equation that is used
as a simple model of phase transition.

Problem 1.
Show that the function
1
v̄(ξ) = (Ex0.0d)
1 + exp(ξ)
satisfies both of the equations
v̄ 0 (ξ) = −v̄(ξ)(1 − v̄(ξ)) (Ex0.0e)
and
v̄ 00 (ξ) = 2v̄(ξ)(1 − v̄(ξ))(1/2 − v̄(ξ)), (Ex0.0f)
and also satisfies
v̄(ξ) → 1 as ξ → −∞,
v̄(ξ) → 0 as ξ → +∞, (Ex0.0g)
v̄(0) = 1/2.

Hence this is a steady state solution for the case α = 1/2.

1
AMATH 586 Hai Zhu Student ID: 1323888

Solution

1
v̄ 0 (ξ) = − eξ
(1 + eξ )2
= −v̄(ξ)(1 − v̄(ξ))

2e2ξ eξ
v̄ 00 (ξ) = −
(1 + eξ )3 (1 + eξ )2

1 2 eξ 1 2 eξ eξ
2v̄(ξ)(1 − v̄(ξ))( − v̄(ξ)) = − +
2 1 + eξ 1 + eξ 2 1 + eξ 1 + eξ 1 + eξ
2e2ξ eξ
= −
(1 + eξ )3 (1 + eξ )2

Thus v̄ 00 (ξ) = 2v̄(ξ)(1 − v̄(ξ))( 21 − v̄(ξ)). Also v̄ satisfies

1 1
v̄(ξ) = → = 1 as ξ → −∞,
1+e ξ 1 + e−∞
1 1
v̄(ξ) = → = 0 as ξ → ∞,
1 + eξ 1 + e∞
1 1
v̄(0) = 0
=
1+e 2
From the above computation, we know v̄ is the solution of

v̄ 00 (ξ) = −2v̄(ξ)(1/2 − v̄(ξ))(v̄(ξ) − 1)

√
Here, ξ = x/δ, δ = 2κ. Thus this is a steady solution for the case α = 1/2.

2
AMATH 586 Hai Zhu Student ID: 1323888

If α 6= 1/2 then the effect of the reaction term is not symmetric. If 0 < α < 1/2 then some
values of v less than 1/2 are driven towards v = 1 by the reaction term. When coupled with
the symmetric diffusion this leads to a traveling wave propagating with some velocity c that
is positive if α < 1/2 or negative if α > 1/2. The traveling wave profile is given by the same
function v̄(x) that satisfies the boundary value problem, and has the form
v(x, t) = v̄((x − ct)/δ), (Ex0.0h)
where the speed c is given by r  
2κ 1
c= −α . (Ex0.0i)
 2

Problem 2.
(a) For any 0 < α < 1, show that v(x, t) = v̄((x − ct)/δ) is a traveling wave solution to
(Ex0.0c) provided that c satisfies (Ex0.0i).
(b) Suppose we define the width of the transition zone (wave front) in a traveling wave to be
the distance in x over which v falls from 0.99 to 0.01. Show that the width of wave front
is roughly 9δ. This can be used to choose a suitable value of h. For example, choosing
h ≈ δ would give roughly 9 grid points in the wave front, which is probably about the
minimum needed to resolve it well numerically. HINT: WLOG, we can pick t = 0 to
answer this since v(x, t) is a travelling wave. So look at v(x, 0) = 1 x .
1+e δ

Solution
(a) Let’s assume we can find a traveling wave solution of the form v̄((x − ct)/δ). Then
substitute into the reaction-diffusion equation, and see if we can find v̄, and whether it
coincides with the solution we get for symmetric case (α = 1/2).
c κ 1
− v̄ 0 = 2 v̄ 00 + v̄(α − v̄)(v̄ − 1) (Ex0.0j)
δ δ 
q √
where c = 2κ  (1/2 − α), δ = 2κ. Actually, we only need c to be some constant times
q
2κ
 , so that we can get rid of other parameters. With a good choice of the constant,
we can also get rid of the effect of α
1 1
− ( − α)v̄ 0 = v̄ 00 + v̄(α − v̄)(v̄ − 1)
2 2
⇒2α(v̄ − v̄(v̄ − 1)) + (v̄ 00 + v̄ 0 − 2v̄ 2 (v̄ − 1)) = 0
0

We already know the solution to v̄ 0 − v̄(v̄ − 1) = 0, which is v̄(ξ) = 1+ξ 1

. Let’s check
00 0 2 0
whether v̄ + v̄ − 2v̄ (v̄ − 1) is also 0 by taking derivative of v̄ − v̄(v̄ − 1) = 0.
0 =(v̄ 0 − v̄(v̄ − 1))0 = v̄ 00 − v̄ 0 (v̄ − 1) − v̄v̄ 0
=v̄ 00 + v̄ 0 − 2v̄v̄ 0 v̄(v̄ − 1) = v̄ 00 + v̄ 0 − 2v̄ 2 (v̄ − 1)

3
AMATH 586 Hai Zhu Student ID: 1323888

Thus v(x, t) = v̄((x − ct)/δ) is a traveling wave solution to (Ex0.0c) provided that c
satisfies (Ex0.0i).

(b) Since v(x, t) is a traveling wave, we can take a look at the transition zone at t = 0, which
is v(x, 0) = 1 x .
1+e δ
Solve v(x, 0) = 0.99, and v(x, 0) = 0.01 gives x1 /δ ≈ −4.595, and x2 /δ ≈ 4.595. Thus
the transition zone is roughly 9δ.

4
AMATH 586 Hai Zhu Student ID: 1323888

Numerical solution The m-file AllenCahn_1d.m solves the equation (Ex0.0c) on the finite
domain a ≤ x ≤ b with no-flux (Neumann) boundary conditions

vx (a, t) = 0, vx (b, t) = 0. (Ex0.0k)

Traveling wave solutions will look similar to those observed for the Cauchy problem until the
wave gets near the boundary.

Problem 3. Download the m-file AllenCahn_1d.m and try it out. You should observe that
the initial discontinuity smooths out into a traveling wave. In parts (a) and (b) make sure
 = .01 is what is set in the code.

(a) Try out both the fzero and newton approaches and observe the difference in efficiency.

(b) Modify the m-file so that it estimates the speed of the traveling wave using the following
observation: For a traveling wave,
Z b
d
v(x, t) dx = c (Ex0.0l)
dt a

as long as the wave front is far from the boundaries. Hence the speed at time tn can be
estimated from numerical results as
 
1  X n X
c≈ h Vj − h Vjn−N  . (Ex0.0m)
Nk
j j

which computes the average speed over the past N time steps (it’s better to take N = 10
or so rather than computing from just one time step so that it behaves more smoothly).
Check that the observed speed approaches a value that agrees with the expected speed,
with an accuracy that is O(k) as the grid is refined. (Refine by the same factor in space
and time.)

(c) Using m = 99 and k = 0.005, change  to  = 0.001. Note that the numerical solution
converges to a traveling wave but that the speed is far from correct, even though the
wave front should be well resolved based on the result of Problem 2(b).
Reduce k with h fixed and note that the error in the speed can be greatly reduced on
this grid by choosing k appropriately.

(d) Try the values m = 99, k = 0.02,  = 0.0001 and explain what you observe. Hint:
What equation is being solved by Newton’s method and how does this behave when k/
is large?

(e) Based on what you observed in parts (b) and (c), explain why we might want to choose
k = min(h, ) to achieve resonable results as we refine the grid.

5
AMATH 586 Hai Zhu Student ID: 1323888

(f) Since we need to take k ≤  to get reasonable solutions to the Allen-Cahn equation,
the reaction ODE is no longer considered to be stiff. This suggests that we could use
an explicit method in this step. Replace backward Euler with the 2-stage Runge-Kutta
method (5.30) and verify that better results are achieved than with the backward Eu-
ler method when stable values of k are used. In fact you should observe second-order
accuracy in the wave speed in this case.

(g) This m-file allows the use of either Crank-Nicolson or TR-BDF2. Since the initial data
is discontinuous, the results of Exercise 9.3 on HW4 suggest that TR-BDF2 may be a
better choice. Find parameter values that illustrate this.

Solution

(a) The following is table of m vs. time when implement fzero

Methods being used: trbdf2 fzero

m CPU time ratio observed order

24.00000 2.13209e-01 NaN NaN
49.00000 7.69770e-01 0.27698 1.79866
99.00000 2.78018e+00 0.27688 1.82594
199.00000 1.05418e+01 0.26373 1.90900

Least squares fit gives T(m) = 0.000600638 * m^1.84221

The following is table of m vs. time when implement newton

Methods being used: trbdf2 newton

m CPU time ratio observed order

24.00000 6.06405e-03 NaN NaN
49.00000 1.59846e-03 3.79369 -1.86803
99.00000 2.39813e-03 0.66654 0.57678
199.00000 4.16304e-03 0.57605 0.78998

Least squares fit gives T(m) = 0.00490876 * m^-0.105655

(b) In order to estimate the speed at time tn , store the results of time tn , and time tn − 10 ∗ k.
Then compute the average speed. The following result is based on ratio h = 5k, tn =
tf inal = 0.3, (speed is given by c in problem 2):

Methods being used: trbdf2 newton

k error ratio observed order

6
AMATH 586 Hai Zhu Student ID: 1323888

0.00500 2.10616e-02 NaN NaN

0.00250 1.08609e-02 1.93921 0.95547
0.00125 5.51343e-03 1.96991 0.97813
0.00063 2.77452e-03 1.98717 0.99071
0.00031 1.38828e-03 1.99852 0.99893

Least squares fit gives E(k) = 3.86159 * k^0.981532

The result by implementing fzero is nearly the same, except that the time is much
longer.

(c) The following is the result by using different value of k:

Methods being used: trbdf2 newton

k error ratio observed order
0.00500 1.67892e+00 NaN NaN
0.00250 3.98168e-01 4.21660 2.07608
0.00125 2.78342e-02 14.30501 3.83845
0.00063 8.97590e-02 0.31010 -1.68920
0.00031 1.38036e-01 0.65026 -0.62092
0.00016 1.53969e-01 0.89652 -0.15760
0.00008 1.63441e-01 0.94205 -0.08613
0.00004 1.68705e-01 0.96879 -0.04574

The following is the log-log plot

Figure 1: error of speed estimate for different k

7
AMATH 586 Hai Zhu Student ID: 1323888

As we can tell from the plot, the error in the speed can be greatly reduced on this grid
by choosing k appropriately.
(d) By choosing  far smaller than k, this numerical method can not give a smooth result of
traveling wave solution. Beyond that, we set α to be 0.3, which make the solution wave
travels to the right. Yet the numerical solution gives a high oscillating result traveling
to the left.
The following is the plot of the final plot:

Figure 2: newton t=0.3, k=0.02, =0.0001

To explain this, let’s first take a look at the behavior when we try to solve backward
Euler.

1 Fvnp = vn + k/epsilon * g(vnp) − vnp;

2 Fprime = k/epsilon * gprime(vnp) − 1;
3 dv = Fvnp ./ Fprime;

If we take a look at the correction vector during each time step, we can find that when
k/ is big, the correction vector can be much bigger than v itself, suppose g(v) is not 0.
This can be verified, if we produce the max norm of the correction vector. We may still
wonder, if newton method doesn’t converge in 10 steps, fzero will be used as a backup
method. So that the solution may still behave well. It turns out for these parameters,
the results are still not close to the real solution even if we use fzero directly.
For this numerical method, let’s consider what it really solves, if we assume the diffusion
and advection part are all solved precisely. For the diffusion part, this numerical method
will give us a solution that will smooth out as time goes by. And for the advection part,
it tries to absorb all the values to the nearest root of the polynomial P . When  is small,
the advection part can be dominant, and the solution will be like

8
AMATH 586 Hai Zhu Student ID: 1323888

Figure 3: fzero t=0.3, k=0.02, =0.0001

If we are implementing newton method, then when we try to get the root for Backward
Euler, it is in fact a combination of newton and fzero. So it hard to tell what it will be
like for a long period of time. But for short time, we can tell from the above two plots
that the results from newton and fezro are similar.
(e) Based on (c), if we choose k much bigger than , then the speed of this traveling wave
can not be well estimated. Based on (b), we usually want to let ratio of k/h be some
constant. Thus it might be a good choice if we let k = min h, .
(f) In order to use 2-stage Runge-Kuttta method, we may modify the code directly by adding
another method after newton and fzero as a replacement of Backward Euler. Or we can
implement 2-stage Runge-kutta first, and then solve the diffusion equation so that the
solution would still satisfy the boundary condition.
The part added is as bellow:

1 case 'rk2'
2 vstar = vn + 1/2*k/epsilon*g(vn);
3 v = vn + k/epsilon*g(vstar);

And the following is the error of estimated speed compared to real speed (ratio h = 5k):

k error ratio observed order

0.00500 4.56999e-03 NaN NaN
0.00250 1.15645e-03 3.95174 1.98249
0.00125 2.97968e-04 3.88112 1.95647
0.00063 8.20948e-05 3.62957 1.85980
0.00031 2.77777e-05 2.95542 1.56336
Least squares fit gives E(k) = 78.2567 * k^1.85405

9
AMATH 586 Hai Zhu Student ID: 1323888

So by changing to 2-stage Runge-Kutta method, we can get second order accuracy when
estimate wave speed. The result after exchanging the order of diffusion and reaction is
still pretty much the same, which would not be surprise, because the wave is still far
from the boundary.
And the following is the log-log plot:

Figure 4: error vs. k

(g) The following is just an example that shows TR-BDF2 is a better choice than CN. Here
I use tfinal= 0.1, k = 0.05, h = 0.03.( = 0.01, κ = 0.3, α = 0.3) So what we will see is
the result after 2 time steps:

Figure 5: result of Crank-Nicolson

10
AMATH 586 Hai Zhu Student ID: 1323888

Figure 6: result of TR-BDF2

11
AMATH 586 Hai Zhu Student ID: 1323888

Numerical methods for the diffusion equation in 2d. The Allen-Cahn equation in
two space dimensions involves two-dimensional diffusion. Before tackling the reaction-diffusion
equation, first consider the diffusion equation. The m-file heat_2d_noflux.m solves the heat
equation in the rectangular domain [ax , bx ] × [ay , by ] in two space dimensions with no-flux
(insulated) boundary conditions at each boundary,

vt (x, y, t) = κ(vxx (x, y, t) + vyy (x, y, t)),

(Ex0.0n)
vx (ax , y, t) = vx (bx , y, t) = vy (x, ay , t) = vy (x, by , t) = 0.

An mx × my grid is used with (mx + 2)(my + 2) unknowns.

Two different methods are implemented:

• If diffmethod = ’CN_unsplit’ then the unsplit Crank-Nicolson method is used, with

the 5-point Laplacian. In this case a linear system of dimension (mx + 2)(my + 2) must
be solved in every time step.

• If diffmethod = ’CN_LOD’ then the Locally One Dimensional method is used based on
dimensional splitting. In this case only tridiagonal systems of dimension mx + 2 and
my + 2 must be solved in each time step.

Problem 4.

(a) In this code the Neumann boundary conditions are imposed using an idea similar to the
“second approach” on page 31 of the text, applied to the method of lines equation at
each grid point that lies on a boundary. Explain why this leads to the modification of
the matrix elements Tx(1,2), Tx(mx+2,mx+1), Ty(1,2), and Ty(my+2,my+1) seen in
the code for the unsplit Crank-Nicolson method.

(b) What boundary conditions are being imposed on U ∗ in the implementation of the LOD
method?

(c) For both the unsplit and LOD versions, test the code on an m × m grid in a square
domain, with different values of m = 24, 49, 99, 199.
Note that the code is written to compute the error and CPU time. Create log-log plots
of both the error for different m and the CPU time required for different m. You should
find that both methods are second order accurate. Does the splitting used in the LOD
method degrade the accuracy very much?
You should find that the unsplit method, which uses Gaussian elimination to solve the
large sparse linear system, is quite slow compared to the LOD method. Suppose that for
large m the CPU time grows like Cmq for some constants C and q. From your computed
results estimate C and q for each method. Approximately how long would it take to
solve this problem with m = 1999 with each method?

(d) Modify the heat_2d_noflux.m code to also implement the TR-BDF2 method in LOD
form.

12
AMATH 586 Hai Zhu Student ID: 1323888

Solution

(a) For the unsplit method, the code is dealing with the following vector:
 
v1,1 v1,2 · · · v1,my +2 · · · vmx +2,1 vmx +2,2 · · · vmx +2,my +2

This can be verified by

1 A = kappa * (kron(Ix,Ty) + kron(Tx,Iy));

For CN method, the discretized form is

Uijn+1 − Uijn κ 2 n
= (D U + Dy2 Uijn + Dx2 Uijn+1 + Dy2 Uijn+1 )
k 2 x ij
Let’s take a look at what happens at the boundary,
1
(v1,j − v−1,j ) = 0
2dx
1
(v1,j − 2v0,j + v−1,j ) = Dx2 U0j
2dx2
Eliminate v−1,j , we can get

1
Dx2 U0j = (−2v0,j + 2v1,j )
2dx2
Thus in the code we need to make modification of Tx(1,2). And once we define the
(mx + 2) × (my + 2) matrix, kron(Tx,Iy) will  lead to the modification
 of the my + 3
to 2(my + 2) columns, which
 is exactly where v2,1 v2,2 · · · v2,my +2 would multiply.
v2,1 v2,2 · · · v2,my +2 are the points next to x = 0 boundary, that need our concern.
The same thing happens to Tx(mx+2,mx+1), Ty(1,2), and Ty(my+2,my+1).

(b) The boundary conditions are being imposed on U ∗ is

vx∗ (ax , y, t) = vx∗ (bx , y, t) = vy∗ (x, ay , t) = vy∗ (x, by , t) = 0.

Let’s take a look at code

1 % x−sweeps.
2 % Note that A2x is applied to all columns of the array v,
3 % and then the system A1x * v = rhs is solved for each column.
4 rhs = A2x*v;
5 v = A1x\rhs;
6
7 % y−sweeps.
8 % in second step, apply on rows:
9 rhs = A2y*v';
10 v = (A1y\rhs)';

13
AMATH 586 Hai Zhu Student ID: 1323888

The x-sweeps says vx∗ (ax , y, t) = vx∗ (bx , y, t) = 0. Because when we deal with v1,j
∗ and
∗
vmx +2,j , matrix A1x has been done the same kind of modification as in (a).
The y-sweeps says vy∗ (x, ay , t) = vy∗ (x, by , t) = 0. Because when we deal with vj,1
∗ and
∗
vj,m , matrix A2y has been done the same kind of modification as in (a).
y +2

This may lead us to wonder that after these two decoupled step, vx∗ (ax , y, t) = vx∗ (bx , y, t) =
0 is not satisfied. Because every value on x-slice has been changed, there is no guarantee
that the derivative with respect to x is still 0 on the boundary. That depends on how we
look at this problem, the value of v0,j∗ is unknown, where v ∗ to v ∗
1,j mx +2,j are grid points.
∗ ∗
We are allowed to still set v0,j = v2,j , which would make it possible to compute next time
step.
In fact, we are assuming that the value of the outside circle is changing according to
pattern of what this LOD method behaves.

(c) The following are the results of unsplit method.

m error ratio observed order

24.00000 3.92507e-03 NaN NaN
49.00000 9.80618e-04 4.00265 -1.94315
99.00000 2.45113e-04 4.00068 -1.97137
199.00000 6.12756e-05 4.00017 -1.98563

Least squares fit gives E(m) = 2.05193 * m^-1.96699

m CPU time ratio observed order

24.00000 6.25219e-02 NaN NaN
49.00000 6.10860e-01 0.10235 3.19341
99.00000 5.72901e+00 0.10663 3.18276
199.00000 4.66254e+01 0.12287 3.00293

Least squares fit gives T(m) = 3.04832e-06 * m^3.13276

Figure 7: unsplit method error vs. m Figure 8: unsplit CPU time vs. m

14
AMATH 586 Hai Zhu Student ID: 1323888

The following are the results of LOD

m error ratio observed order

24.00000 3.93780e-03 NaN NaN
49.00000 9.83832e-04 4.00251 -1.94310
99.00000 2.45919e-04 4.00064 -1.97136
199.00000 6.14772e-05 4.00016 -1.98563

Least squares fit gives E(m) = 2.05847 * m^-1.96697

m CPU time ratio observed order

24.00000 2.94121e-03 NaN NaN
49.00000 1.68664e-02 0.17438 2.44688
99.00000 1.19709e-01 0.14089 2.78650
199.00000 9.42140e-01 0.12706 2.95493

Least squares fit gives T(m) = 4.52111e-07 * m^2.7328

Figure 9: LOD method error vs. m Figure 10: LOD CPU time vs. m

As we can see, the LOD method is still second order accurate.

By using the result polynomial T (m) = 3.04832e − 06 ∗ m3 .13276 for unsplit,

T (1999) = 3.04832e − 06 ∗ 19993.13276 = 6.6791e + 04

The time estimate for splitting method is

T (1999) = 4.52111e − 07 ∗ m2.7328 = 4.7382e + 03

(d) The following is the modified part for implementing TR-BDF2 when we try to set up
some necessary matrix:

15
AMATH 586 Hai Zhu Student ID: 1323888

1 switch diffmethod
2
3 case 'TRBDF2 LOD'
4
5 ex = ones(mx+2,1);
6 Tx = spdiags([ex −2*ex ex],[−1 0 1],mx+2,mx+2);
7 Tx(1,2) = 2;
8 Tx(mx+2,mx+1) = 2;
9 Tx = kappa/dxˆ2 * Tx;
10 % matrix for step1 CN
11 A1x = speye((mx+2)) − k/4 * Tx;
12 A2x = speye((mx+2)) + k/4 * Tx;
13 % matrix for step2 BDF
14 A3x = speye(mx+2) − k/3 * Tx;
15
16 ey = ones(my+2,1);
17 Ty = spdiags([ey −2*ey ey],[−1 0 1],my+2,my+2);
18 Ty(1,2) = 2;
19 Ty(my+2,my+1) = 2;
20 Ty = kappa/dyˆ2 * Ty;
21 % matrix for step1 CN
22 A1y = speye((my+2)) − k/4 * Ty;
23 A2y = speye((my+2)) + k/4 * Ty;
24 % matrix for step2 BDF
25 A3y = speye(my+2) − k/3 * Ty;

The following is the modified part for implementing TR-BDF2 when we compute during
one time step:

1 switch diffmethod
2
3 case 'TRBDF2 LOD'
4
5 % Use TR−BDF2 in LOD form:
6
7 % x−sweeps.
8 % Note that A2x is applied to all columns of the array v,
9 % and then the system A1x * v = rhs is solved for each column.
10 % TR step:
11 rhs = A2x*vn;
12 vnph = A1x\rhs;
13 % BDF2 step:
14 rhs = (4*vnph − vn) / 3;
15 vn = A3x\rhs;
16
17 % y−sweeps.
18 % in second step, apply on rows:
19 % TR step:
20 rhs = A2y*vn';
21 vnph = (A1y\rhs)';
22 % BDF2 step:
23 rhs = (4*vnph − vn)' / 3;
24 vn = (A3y\rhs)';

16
AMATH 586 Hai Zhu Student ID: 1323888

The following is the results for different values of m:

m error ratio observed order

24.00000 3.95476e-03 NaN NaN
49.00000 9.88159e-04 4.00215 -1.94298
99.00000 2.47006e-04 4.00055 -1.97132
199.00000 6.17494e-05 4.00013 -1.98562

Least squares fit gives E(m) = 2.06703 * m^-1.96691

m CPU time ratio observed order

24.00000 6.21131e-03 NaN NaN
49.00000 3.31999e-02 0.18709 2.34835
99.00000 2.50243e-01 0.13267 2.87201
199.00000 1.93287e+00 0.12947 2.92806

Least squares fit gives T(m) = 9.4445e-07 * m^2.72891

Figure 11: TR-BDF2 in LOD form er- Figure 12: TR-BDF2 in LOD form
ror vs. m CPU time vs. m

17
AMATH 586 Hai Zhu Student ID: 1323888

The two-dimensional Allen-Cahn equation. In two space dimensions the Allen-Cahn

equation takes the form
1
vt = κ(vxx + vyy ) + g(v), (Ex0.0o)

where g(v) is again given by (Ex0.0b).
If α = 1/2 then an interface between two materials is stationary only if it is flat. It can be
shown that a curved interface moves normal to the interface at a velocity that is proportional
to the curvature at each point.
Please see the Sample Results on the webpage,
[Link]
to see the behavior when the initial data is u = 1 inside a rectangle and u = 0 outside. The
square corners become rounded and the shape contracts towards a circular shape before finally
disappearing. In fact the same things happens for any initial shape — it smooths out and
ultimately contracts as a small circular shape towards a vanishing point. This behavior is
similar to what’s seen as an ice cube melts.

Problem 5. Combine the techniques from the m-files used above to create a program that
solves the two-dimensional Allen-Cahn equation in a general rectangular domain with no-flux
boundary conditions.
For diffusion use the LOD form of TR-BDF2. For the ODE solver use the 2-stage Runge-
Kutta method (5.30).
Test your code by solving the following problems with κ = 0.1, α = 0.5,  = 0.001 on a
100 × 100 grid in −1 ≤ x ≤ 1, − 1 ≤ y ≤ 1:

(a) The example shown on the webpage above under Sample Results, which can be specified
by

v = zeros(size(X));
v(X>-0.8 & X<0.8 & Y>-0.3 & Y<0.3) = 1.0;

(b) Using initial data

v(x, y, 0) = 0.2 1 + (1 + cos(2πx))(1 + cos(4πy)) . (Ex0.0p)

With this data you should see a periodic array of spots.

(c) With data


 1 if − 0.5 < x < 0, − 0.5 < y < 0,
v(x, y, 0) = 1 if − 0.1 < x < 0.8, − 0.1 < y < 0.8, (Ex0.0q)
0 otherwise


Print out plots of a few times for each problem to show that your code is working.

18
AMATH 586 Hai Zhu Student ID: 1323888

Solution
The LOD form of TR-BDF2 is the same as in problem 4(d). The ODE solver is also
the same as in problem 3(f). The part that contains the initial conditions is as follows:

1 % set up initial data

2 switch initialdata
3
4 case 'rectangular'
5 v = zeros(size(X));
6 v(X>−0.8 & X<0.8 & Y>−0.3 & Y<0.3) = 1.0;
7
8 case 'periodic'
9 v = 0.2*(1 + (1 + cos(2*pi*X))).*(1 + cos(4*pi*Y));
10
11 case '2rectangular'
12 v = zeros(size(X));
13 v(X>−0.5 & X<0.0 & Y>−0.5 & Y<0.0) = 1.0;
14 v(X>−0.1 & X<0.8 & Y>−0.1 & Y<0.8) = 1.0;
15 end % switch

(a) The following are the results for rectangular initial condition:

Figure 13: solution of t = 0 Figure 14: solution of t = 0

19
AMATH 586 Hai Zhu Student ID: 1323888

Figure 15: solution of t = 1 Figure 16: solution of t = 1

Figure 17: solution of t = 8 Figure 18: solution of t = 8

(b) The following are the results for periodic initial condition:

Figure 19: solution of t = 0 Figure 20: solution of t = 0

20
AMATH 586 Hai Zhu Student ID: 1323888

Figure 21: solution of t = 1 Figure 22: solution of t = 1

Figure 23: solution of t = 2 Figure 24: solution of t = 2

(c) The following are results for two rectangular:

Figure 25: solution of t = 0 Figure 26: solution of t = 0

21
AMATH 586 Hai Zhu Student ID: 1323888

Figure 27: solution of t = 1 Figure 28: solution of t = 1

Figure 29: solution of t = 3 Figure 30: solution of t = 3

22
AMATH 586 Hai Zhu Student ID: 1323888

Extra Credit Problems

1
v 0 (t) = (g(v(t)) − w(t) + Ia ),
 (Ex0.0r)
w0 (t) = βv(t) − γw(t).

Problem 6. (5 points)

(a) Write a code using ode15s to solve (Ex0.0r) and reproduce the figures shown above as a
test that it is working.

(b) Modify your code to also produce a phase plane plot (i.e., a plot of w(t) vs. v(t) as they
vary in time) similar to the plots on
[Link]
for example. Do this for each set of parameters used above, i.e. for

Ia = 0, v0 = 0.29,
Ia = 0, v0 = 0.31,
Ia = 0.2, v0 = 0.

(c) Experiment with varying  and comment on what you observe, both in terms of how
the solution behaves (as a function of t and in the phase plane) and in terms of the
computational method.

Solution

(a) The following is the main part of my code:

1 [t,v] = ode15s(@myfun, [0 T], [v0(1), v0(2)]);

2
3 function vp = myfun( t, v0)
4 v = v0(1);
5 w = v0(2);
6 vp1 = (v*(alpha − v)*(v − 1) − w + I)/epsilon;
7 vp2 = beta*v − gamma*w;
8 vp = [vp1; vp2];
9 end

The following are the figures:

23
AMATH 586 Hai Zhu Student ID: 1323888

Figure 31: solution of v0 = 0.29, Iα = Figure 32: solution of v0 = 0.31, Iα =

0 0

Figure 33: solution of v0 = 0, Iα = 0.2

(b) The following are the phase plot for different parameters:

Figure 34: solution of v0 = 0.29, Iα = 0

24
AMATH 586 Hai Zhu Student ID: 1323888

Figure 35: solution of v0 = 0.31, Iα = 0

Figure 36: solution of v0 = 0, Iα = 0.2

(c) By changing values of , we can observe the solution barely change even if we change
the  enormously. And when it comes to phase, I think it behaves the same way, which
means the ratio of w(t) vs. v(t) doesn’t change to much. The following is the computing
time for different :

time spend with epsilon= 0.000001 is 5.7197e-02

time spend with epsilon= 0.00001 is 5.5911e-02
time spend with epsilon= 0.0001 is 4.4839e-02

25
AMATH 586 Hai Zhu Student ID: 1323888

time spend with epsilon= 0.001 is 4.9097e-02

The following is the figure for different , it looks like we can not distinguish one from
another.

Figure 37: solution for different  v0 = 0.31, Iα = 0

Figure 38: phase plot for different  v0 = 0.31, Iα = 0

26
AMATH 586 Hai Zhu Student ID: 1323888

Problem 7. (5 points)

Write a Matlab m-file to solve the one-dimensional FitzHugh-Nagumo equations (??). You
can base your code on AllenCahn_1d.m. You will have to keep track of w as well as v and
modify the reaction terms for the FitzHugh-Nagumo reactions. Again use the backward Euler
method for these terms. This will now require solving a system of 2 equations for Vjn+1 and
Wjn+1 at every grid point. Note, however, that the equation for w is linear and so it is possible
to express Wjn+1 as a function of Vjn+1 . Use this to reduce the problem to a scalar cubic
equation to be solved for Vjn+1 and use Newton’s method to solve this, similar to the approach
used in AllenCahn_1d.m. (You might want to first try using fzero at every point since this
is simpler, but you will find it runs very slowly.) Also implement the 2-stage explicit Runge-
Kutta method (5.30) and comment on whether this is an efficient approach relative to the
implicit method.
Test your code with the following two tests (feel free to experiment with more):

(a) Use m = 200, k = 0.002 on 0 ≤ x ≤ 6 with parameter values

α = 0.3, β = 1, γ = 1, κ = 0.2,  = 0.001 (Ex0.0s)

with no applied current (Ia = 0) and initial data


1 if x < 1
v(x, 0) = , w(x, 0) = 0.
0 if x > 1

You should see a traveling wave develop and propagate.

(b) Same situation as in part (a) but with initial data v(x, 0) = w(x, 0) = 0 and an applied
current Ia (x) = 0.8 exp(−5x2 ). This models a situation in which a stimulated nerve cell
sends out a train of pulses. Run the computation out to time t = 3 or so and you should
see several pulses generated.

Solution
The modified part of code is as follows:

1 %I = @(x) 0.8*exp(−5*x.ˆ2);
2 I = @(x) 0;
3
4 switch odemethod
5 case 'rk2'
6 vstar = vn + 1/2*k/epsilon*(g(vn)...
7 ...−(k*beta)*vnp/(1+k*gamma)−wn/(1+k*gamma)+I(x));
8 v = vn + k/epsilon*(g(vstar)−...
9 ...(k*beta)*vstar/(1+k*gamma)−wn/(1+k*gamma)+I(x));
10 w = (k*beta*v+wn)/(1+k*gamma);
11 case 'newton'
12 for inewton=1:10
13 Fvnp = vn + k/epsilon * (g(vnp)...
14 ...−(k*beta)*vnp/(1+k*gamma)−wn/(1+k*gamma)+I(x)) − vnp;
15 Fprime = k/epsilon * (gprime(vnp)−k*beta/(1+k*gamma)) − 1;

27
AMATH 586 Hai Zhu Student ID: 1323888

16 dv = Fvnp ./ Fprime;
17 vnp = vnp − dv;
18 if maxF < 1e−6
19 % convergence at all grid points
20 break
21 end
22 end
23
24 v = vnp;
25 ...
26 w = (k*beta*v+wn)/(1+k*gamma);

We don’t need to change the diffusion part, because w doesn’t have diffusion. And the
time difference between newton and 2-stage Runge-Kutta are as follows:

1 Methods being used: trbdf2 newton

2
3 time spend on diffusion = 8.3729e−03
4 time spend on ODEs = 2.7679e−02
5
6 Methods being used: trbdf2 rk2
7
8 time spend on diffusion = 7.1074e−03
9 time spend on ODEs = 1.4009e−02

(a) The following are the figures of this traveling wave

Figure 39: solution of t = 0 Figure 40: solution of t = 0.1

28
AMATH 586 Hai Zhu Student ID: 1323888

Figure 41: solution of t = 0.15 Figure 42: solution of t = 1

(b) The following are the figures of this traveling wave of different parameter:

Figure 43: solution of t = 0.15 Figure 44: solution of t = 1

Figure 45: solution of t = 1.5 Figure 46: solution of t = 3.3

29