# 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, 2016 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, 2016 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, 2016 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, 2016 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, 2016 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, 2016 at 13:13
• @AlexeiBoulbitch spinodal decomposition needs random initial value t=0 . Jul 31, 2019 at 21:55