Solving the nonlinear Cahn-Hilliard equation

I am trying to solve the Cahn-Hilliard equation with a "random" initial condition and boundary conditions as described below:

$$h_t = \nabla^2\left(-\gamma \nabla^2 h + h^3 - h\right)$$ Boundary conditions:

$$h^{(1,0,0)}(0,y,t)=0,h^{(1,0,0)}(L,y,t)=0$$ $$h^{(3,0,0)}(0,y,t)=0,h^{(3,0,0)}(L,y,t)=0$$ $$h^{(0,1,0)}(x,0,t)=0,h^{(0,1,0)}(x,L,t)=0$$ $$h^{(0,3,0)}(x,0,t)=0,h^{(0,3,0)}(x,L,t)=0$$

For now, I am using a "sine wave" initial condition (yes, the initial and boundary conditions are not compatible). My Mathematica code for this is as follows:

Needs["DifferentialEquationsInterpolatingFunctionAnatomy"];
Clear[Eq5, Complete]
Eq5[h_, {γ_}] := \!$$\*SubscriptBox[\(∂$$, $$t$$]h\) -
Laplacian[-γ Laplacian[h, {x, y}] + h^3 - h, {x, y}] == 0;
SetCoordinates[Cartesian[x, y, z]];
Complete[γ_] := Eq5[h[x, y, t], {γ}];

L = 1.; TMax = 1.;
bsf = Interpolation@
Flatten[Table[{{x, y}, 1 + .05*RandomReal[{-1, 1}]}, {x, 0,
L + 1}, {y, 0, L + 1}], 1];
Off[NDSolve::mxsst];
Off[NDSolve::ibcinc];
hSol = h /. NDSolve[{
Complete[0.001],
h[x, y, 0] == 1 + 0.5 Sin[2 π x/L] Cos[2 π y/L],
(*h[x,y,0]\[Equal]bsf[x,y],*)
Derivative[1, 0, 0][h][0, y, t] == 0,
Derivative[1, 0, 0][h][L, y, t] == 0,
Derivative[3, 0, 0][h][0, y, t] == 0,
Derivative[3, 0, 0][h][L, y, t] == 0,

Derivative[0, 1, 0][h][x, 0, t] == 0,
Derivative[0, 1, 0][h][x, L, t] == 0,
Derivative[0, 3, 0][h][x, 0, t] == 0,
Derivative[0, 3, 0][h][x, L, t] == 0
},
h,
{x, 0, L},
{y, 0, L},
{t, 0, TMax}, Method -> "LSODA"
][]

This takes forever to run (>5 min on a machine with 32 gigs of ram, I quit it). Are there any tuning parameters that I should use for this equation for it to run smoothly without warnings? It is a stiff equation and hence the choice of LSODA (mma would have chosen this automatically anyway?)

I do notice that for negative values of $\gamma$, stiffness is arrived at sooner and code terminates within 1-2 seconds on my computer.

I also have the warning: "Requested order is too high; order has been reduced to {2,2}."

• bsf isn't used actually ? Are you still in v8 and using VectorAnalysis` package? Sep 7 '16 at 3:44
• @xzczd I didn't use bad because that initial condition could have exacerbated the stiffness in the problem. But yes I must use it. Also, I am using (explicitly calling) VectorAnalysis out of force of habit! Sep 7 '16 at 13:00

I tried a slightly different boundary conditions, mainly since the solution with this conditions is easier:

Clear[s, eq, bc1, bc2, ic, Lx1, Ly, \[Eta], x, y, t];
Ly = 1;
Lx = 1;

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

bc1P = \[Eta][t, x, -Ly] == \[Eta][t, x, Ly];
bc2P = \[Eta][t, -Lx, y] == \[Eta][t, Lx, y];
ic = \[Eta][0, x, y] == 0.5*Exp[-4 (x^2 + y^2)];

and the MethodOfLines

mol = {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid",
"DifferenceOrder" -> "Pseudospectral"}};

sol = NDSolve[{eq, bc1P, bc2P, ic}, \[Eta], {t, 0, 1}, {x, -Lx,
Lx}, {y, -Ly, Ly}, Method -> mol] // AbsoluteTiming

It took 703.5 s, but the solution was obtained. Here is its visualization

Grid@Partition[
Table[Plot3D[\[Eta][T, x, y] /. sol[[2, 1, 1]], {x, -Lx,
Lx}, {y, -Ly, Ly}, PlotRange -> All,
AxesLabel -> {"x", "y", "\[Eta]"}, PlotPoints -> 31,
PlotLabel ->
Row[{Style["T=", Italic, 12], Style[T, Italic, 12]}]], {T, {0,
0.005, 0.01, 0.05, 0.07, 0.09, 0.1, 0.2, 0.3, 0.4, 0.5, 1}}], 3]

yielding this: I hope it helps. On the other hand the solution does not look like a spinodal decomposition, at least at the first glance.

Have fun!

• What do you recommend for spinodal decomp to be seen? It shouldn't be the initial conditions, right? Irrespective of the IC, shouldn't the same attractor be attained? Sep 7 '16 at 11:47
• @drN I do not know. I never worked on it myself, but observed a number of experiments. And what I have seen looked differently. Besides, I am not quite sure, but is h in the equation not a concentration? If yes, one needs to account that it must be positive. Further, I know that there are some accepted theories describing at least initial stages of the decomposition. May be one needs to start from them and simulate their results. But very frankly, I do not believe that the homogeneous Cahn-Hilliard really describes any realistic situation. But this subject comes out of the scope of this site Sep 7 '16 at 12:39
• Yes, h could be concentration if we consider the Cahn-Hilliard formulation. However, the CH equation is a special case that can be reduced (under certain conditions) to other evolutionary equations. I found a MATLAB code that does the CH and arrives at spinodal decomposition. Interesting. Sep 7 '16 at 13:02
• @drN Ok, as I said, I never looked into it precisely, and do not know the state of the art. People work on it already several decades, some of them can program. So I think one will for sure find some numerical approaches. Sep 7 '16 at 13:13
• @AlexeiBoulbitch spinodal decomposition needs random initial value t=0 . Jul 31 '19 at 21:55