# Numerical solution of the 2D-spatial nonlinear Allen equation

I would like to solve the 2D-spatial Allen equation in rectangular coordinate, which is a nonlinear reaction-diffusion PDE of the type

$$\partial_{t}u=\epsilon(\partial_{xx}+\partial_{yy})u + u - u^{3},$$

subject the initial condition $$u(x,y,0)=sin(2 \pi x)+0.001cos(16 \pi x)$$, where $$\epsilon$$ is a small but positive constant.

Let's now take the FFT (which include periodic boundary condition) of both sides of the Allen equation to obtain

$$\widehat{\partial_{t}u}_{k} = \epsilon\,(\widehat{\partial^{2}_{xx}u_{k}}+\widehat{+\partial^{2}_{xx}u_{k}})+\widehat{u}_{k}-\widehat{u^{3}}_{k}$$

Where $$\widehat{\partial^{2}_{xx}u_{k}}=(i\,k_{x})^{2}\widehat{u}_{k}$$ and $$\widehat{\partial^{2}_{yy}u_{k}}=(i\,k_{y})^{2}\widehat{u}_{k}$$.

Taken into account that nonlinear term $$FFT(u^{3}) \neq FFT(u)^3$$, so the $$u^{3}$$ must be computed before taking the FFT. Then,

$$\partial_{t}\widehat{u}_{k} = \epsilon\,((ik_{x})^{2}+(ik_{y})^{2})\,\widehat{u}_{k}+\widehat{u}_{k}-\widehat{u^{3}}_{k}$$

In order to solve this numerically we are going to use a combination of implicit (backward Euler) and explicit (forward Euler) methods (Euler is unstable for this equation). Applying this to Allen equation we find that

$$\widehat{u}^{n+1}_{k} = \frac{\widehat{u}\,^{n}_{k}(1+1/h)+\widehat{u}_{k}-\widehat{(u^{n})^{3}}_{k}}{-\epsilon\,(ik_{x})^{2}-\epsilon\,(ik_{y})^{2}+1/h}$$

where $$k_{x}$$ and $$k_{y}$$ is to remind us that we take the FFT in respected directions. Notice that when programming we are going to have to update the nonlinear term $$(u^{3})$$ each time you want to calculate the next timestep $$n+1$$. The reason this is worth mentioning is that for each timestep we are going to have to go from real space to Fourier space to real space, then repeat.

Edit: There is a version of the code in Matlab here.

n = 500;(*Number of discretization points*)
T = 10; (*Time Integration*)
\[Epsilon] = 10^-2; (*diffusivity*)
dt = 10^-2; (*step time*)
kx = I Join[Range[0, n/2 - 1], {0},Range[-n/2 + 1, -1]];(* ik vector in x direction*)
ky = I Join[Range[0, n/2 - 1], {0},Range[-n/2 + 1, -1]];(* ik vector in y direction*)
k2x = kx^2; (*(ik)^2 vector in x direction*)
k2y = ky^2;(*(ik)^2 vector in x direction*)
icu = Table[Sin[2 Pi i] + 0.001 Cos[16 Pi i], {i, n}, {j,n}]; (*initial condition to u*)
ictu = Fourier[icu,FourierParameters -> {1, -1}];(*%FFT for linear*)
icv = Table[Sin[2 Pi i] + 0.001 Cos[16 Pi i], {i,n}, {j,n}]; (*initial condition to u*)
var = Table[{Subscript[u, i, j][t], Subscript[v, i, j][t]}, {i,n}, {j,n}]; (*vector of variables*)
var = Flatten[var]; (*vector of variables*)
eqns = Table[{Subscript[u, i, j][t + 1] == (Subscript[u, i, j][t] (1/dt + 1) - (Subscript[v, i,j][t])^3)/(-\[Epsilon] (k2x[[i, j]] + k2y[[i, j]]) + 1/dt), Subscript[u, i, j] == ictu[[i, j]], Subscript[v, i, j] == ictv[[i, j]]}, {i, n}, {j, n}];
eqn = Flatten[eqns];
sol = RecurrenceTable[eqn,vars, {t, 0, T}];(*Simulate in real and Fourier frequency domain*)


However I am not getting to update the nonlinear term $$(u^{3})$$ each time you want to calculate the next timestep n+1. Can anybody help me?

a = Table[InverseFourier[su[[i]]], {i, Length[su]}];
b = Table[{j, l, Re[a[[i, j, l, 1]]]},{i,Length[a]}, {j,Length[a[[i]]]}, {l, Length[a[[i,j]]]}];
c = Table[Flatten[b[[i]], 1], {i, Length[b]}];
Table[ListPlot3D[c[[i]], PlotRange -> All], {i, Length[c]}]


I'm looking for a solution in stationary state like • Perhaps Method->"FiniteElement"solves your problem too! Nov 25, 2021 at 13:16
• Is it still for didactic purpose? Are you aware that the method is aleady implemented in NDSolve as Pseudospectral method? Nov 25, 2021 at 13:36
• Hi @xzczd, yes, it is for didactic purposes. I implemented it using NDSolve with Pseudospectral
– SAC
Nov 25, 2021 at 14:58
• \[Epsilon] = 0.001; Ly = 1; Lx = 1; eq = D[u[t, x, y], t] == Laplacian[u[t, x, y], {x, y}] + u[t, x, y] - u[t, x, y]^3; bc1P = u[t, x, -Ly] == u[t, x, Ly]; bc2P = u[t, -Lx, y] == u[t, Lx, y]; ic = u[0, x, y] == Sin[2 Pi x] + 0.001 Cos[16 Pi x]; mol = {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "DifferenceOrder" -> "Pseudospectral"}}; sol = NDSolve[{eq, bc1P, bc2P, ic}, u, {t, 0, 1}, {x, -Lx, Lx}, {y, -Ly, Ly}, Method -> mol], but it doesn't seem to work.
– SAC
Nov 25, 2021 at 15:00
• 1. What do you mean by "it doesn't seem to work"? 2. Spatial grid is too coarse. 3. End of time is too small. 4. \[Epsilon] is missing in your equation. Nov 25, 2021 at 15:26

First, I would like to make two notes.

1. This equation has also another very famous name: the time-dependent Ginzburg-Landau equation. As such, it was studied in a huge number of contexts.

2. The parameter eps should be removed from the very beginning by rescaling x->Sqrt(eps) x; y->Sqrt(eps) y. After that, we come to the following equation:

 eq = D[u[t, x, y], t] ==
Laplacian[u[t, x, y], {x, y}] + u[t, x, y] - u[t, x, y]^3;


This equation, as it is written down here describes a relaxation of the dependent variable u to its equilibrium values u=1 or u=-1.

During the numerical solution of this equation one often faces a trap. Under some boundary or initial conditions, the solution relaxes to u=0 instead of u=1 or u=-1. This is related to the peculiarities of the implemented numeric procedure.

This, for example, happens when the boundary conditions u(t,boundary)=0 are used. The same can happen if the initial conditions are symmetrical with respect to the planes u=0 and, say, x=0. This happens in the case of your initial conditions.

There can be several approaches to overcome this trap. The right approach depends on the nature of the problem that you solve. For example, if applicable, let us use an asymmetric initial condition:

ic = u[0, x, y] == Sin[2 Pi (x - 1/3)] + Exp[-(x - 1/3)];


With the same other settings:

Ly = 1;
Lx = 1;
bc1P = u[t, x, -Ly] == u[t, x, Ly];
bc2P = u[t, -Lx, y] == u[t, Lx, y];
mol = {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid",
"DifferenceOrder" -> "Pseudospectral"}};


sol = NDSolve[{eq, bc1P, bc2P, ic},
u, {t, 0, 5}, {x, -Lx, Lx}, {y, -Ly, Ly}, Method -> mol][[1, 1]]


The relaxation one observes as follows:

Animate[Plot3D[
Evaluate[u[t, x, y] /. sol], {x, -Lx, Lx}, {y, -Ly, Ly},
PlotRange -> {0, 3.5}], {t, 0, 0.10}, AnimationRate -> 0.01,
AnimationRepetitions -> 1]


Below is the initial stage of the evolution next comes an intermediate step: and the final step of simulation at t=0.3: Have fun!

• "The parameter eps should be removed from the very beginning by rescaling x->Sqrt(eps) x; y->Sqrt(eps) y. " Well, it's not necessary for obtaining a numeric solution, isn't it? Nov 26, 2021 at 0:21
• @xzczd It is not necessary to be able to solve this simplified equation, you are right. However, first, the possibility to remove it shows that it only influences the solution scale, but not its structure. Second, in practice one usually deals with the equation with coefficients: kappa*D[u[x,y,t],t]==g*D[u[t,x,y],x,x]+a*u[x,y,t]-b*u[x,y,t]^3. In this case, I would recommend rescaling the equation such as removing as many coefficients as possible. Last, a small factor in front of the higher derivative has the potential to make some difficulties. Nov 26, 2021 at 11:13
• @xzczd Continuation: For example, if there is kink-like behavior, the kink width is ~Sqrt[eps], and with small eps it will require a smaller mesh size within the kink. However, in general, I agree with you that rescaling is not absolutely necessary here. Nov 26, 2021 at 11:16
• @AlexeiBoulbitch can this "mol" be applied for fourth order differential equations (e.g.cahn hilliard)? Dec 15, 2021 at 1:33
• @ABCDEMMM "mol" is simply the name of the set of methods. As such one can apply it to any suitable equation. However, as much as I know, so far no possibility to directly solve equations of the fourth-order has been implemented. All cases of solving such an equation that I have seen required to introduce a new dependent variable such that the equation of the fourth-order would transform into a system of two second-order equations. Then this system can be solved using "mol". For the Cahn-Hillard equation, this transform does not work due to nonlinear term. At least, I do not see how to make it. Dec 15, 2021 at 10:07

## NDSolve

First, the NDSolveValue solution which appeared in the comments above to compare with (edit: added MinPoints->35 per @xzczd's suggestion):

ϵ = 0.001;

eq = D[u[t, x, y], t] == ϵ Laplacian[u[t, x, y], {x, y}] + u[t, x, y] - u[t, x, y]^3;
bc1P = u[t, x, -1] == u[t, x, 1];
bc2P =u[t, -1, y] == u[t, 1, y];
ic = u[0, x, y] == Sin[2 Pi x] + 0.001 Cos[16 Pi x];

sol = NDSolveValue[{eq, bc1P, bc2P, ic},
u, {t, 0, 5}, {x, -1, 1}, {y, -1, 1},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"DifferenceOrder" -> "Pseudospectral", MinPoints -> 35}}] The two domains have 'phase-separated' with a 'diffuse' interface, consistent with Allen-Cahn dynamics.

## Spectral Iterative Solver

For the spectral method we first make the 2D k-grid and precompute k^2.

Note: I think your wavevectors included the origin twice. Here I used the ifftshift function from here to shift an equispaced grid.

ifftshift[dat_?ArrayQ, k : (_Integer?Positive | All) : All] :=
Module[{dims = Dimensions[dat]},
RotateRight[dat,
If[k === All, Ceiling[dims/2],
Ceiling[dims[[k]]/2] UnitVector[Length[dims], k]]]]

wavevector1D[n_] := ifftshift[Rest@Subdivide[-n/2, n/2., n]]
wavevectorSquared[{nx_, ny_}] :=
With[{kx = ConstantArray[wavevector1D[nx], ny]\[Transpose],
ky = ConstantArray[wavevector1D[ny], nx]},
kx^2 + ky^2]


We then implement a single iteration of the time-stepping scheme as:

spectralAllen[u_, {ϵ_,  dt_}][k2_] :=
With[{fourierU = Fourier[u, FourierParameters -> {1, -1}],
fourierUCubed  = Fourier[u^3, FourierParameters -> {1, -1}]},
Re@InverseFourier[(fourierU (1 + 1/dt) -
fourierUCubed)/(1/dt + ϵ k2),
FourierParameters -> {1, -1}]
]


define initial conditions:

dt = 0.001;
n = 128;
uvals = Table[Sin[2 \[Pi] x] + 0.001 Cos[16 \[Pi] x],
{y,Rest[Subdivide[-1, 1, n]]}, {x, Rest[Subdivide[-1, 1, n]]}];
k2 = wavevectorSquared[{n, n}];


and iterate:

Monitor[
checkpointedSolution =
Reap[Do[
uvals = spectralAllen[uvals, {ϵ, dt}][k2];
If[Mod[t, 100] == 0, Sow[uvals]]
, {t, 1, 5000}]][[2, 1]];
, ProgressIndicator[t, {1, 5000}]
]


This will keep updating uvals at each iteration and further Reap/Sow some checkpoints.

The two methods seem to agree reasonably well:

makePeriodicRank2[tensor_] := Block[{array = tensor},
array = Join[Map[List, array[[-1]], {0}], array, 1];
array = Join[Map[List, array[[All, -1]], {1}], array, 2]]

Show[
ListPlot3D[makePeriodicRank2[uvals], DataRange -> {{-1, 1}, {-1, 1}},
PlotStyle -> Red],
Plot3D[sol[5, x, y], {x, -1, 1}, {y, -1, 1}, PlotStyle -> Blue]
] • As mentioned in my comment above, the spatial grid is too coarse. Try MinPoints->35 Nov 26, 2021 at 0:14
• Thanks, updated - NDSolveValue still seems to get 'stuck' though and not 'phase-separating' all the way, even with very fine discretization Nov 26, 2021 at 0:43
• Which version are you in? In v12.3.1 NDSolveValue solves the problem instantly and the result is expected: i.stack.imgur.com/vF5kb.png Nov 26, 2021 at 1:27
• @GeorgeVarnavides great job!
– SAC
Nov 26, 2021 at 3:36
• Pretty sure the $\LaTeX$ formula is wrong (extra factor of $\hat{u}_k$). It seems like the OP is using a simple Forward-Backward Euler (1,1,1) IMEX scheme, see e.g. section 2.1 here. My wavevectors were scaled wrong - updated, solutions are consistent now. Nov 26, 2021 at 15:51