How to make a Cahn–Hilliard model GIF

Question

The question derived from @Jason B 's answer.I wanna make it but don't how to do.But we can get some rule about the GIF

1. The change of the total pixel is very small.

imagelist = Import["https://i.stack.imgur.com/y9np9.gif"]

totalpixel = 
  ImageMeasurements[#, "Total"] & /@ imagelist // Total /@ # &;
ListLinePlot[Differences[totalpixel]/Most@totalpixel, Mesh -> All]

2. A frame texture and the next frame roughly is the same.

So how to make this GIF by its first frame?

Answer 1

I am delighted by this problem mostly because I was not aware of the underlying physical model of phase separation (the Cahn–Hilliard equation)!

Anyway, here is an approximation of a somewhat similar behavior to start the discussion:

NestList[
    Sharpen@Threshold[GaussianFilter[#, 2], {"Hard", {"BlackFraction", 0.1}}] &, 
    img, 50
];

ListAnimate[%]

I am looking forward to something better from the image processing experts!

---

Here is an alternative method. It is still based on image processing, but it uses CommonestFilter instead of the combination of operations used before. It is interesting to note that these situations almost always reach a steady state in relatively few iterations (hence the use of FixedPointList rather than NestList here.

SeedRandom[20160216]
img = Image@RandomChoice[{0, 1}, {256, 256}];
list = FixedPointList[Binarize[CommonestFilter[#, 3], 0.5] &, img, 2000];

Answer 2

Taken from the Matlab code here http://www.math.utah.edu/~eyre/computing/matlab-intro/ch.txt with slight modifications:

m = n = 256;
delx = 1/(m - 1);
delx2 = delx^2;
x = Range[0, 1, delx];

delt = 0.00000001;
ntmax = 250;

epsilon = 0.01;
eps2 = epsilon^2;

a = 2;

lam1 = delt/delx2;
lam2 = eps2*lam1/delx2;

leig = Transpose@Table[2 Cos[π Range[0., n - 1]/(n - 1)] - 2, {m}] + 
 Table[2 Cos[π Range[0., m - 1]/(m - 1)] - 2, {n}];

cheig = 1 - (a*lam1*leig) + (lam2*leig^2);

seig = lam1*leig;

u = RandomReal[{-.5, .5}, {n, m}];
hatu = FourierDCT[u];

t = 0.;
Monitor[
 While[t < ntmax*delt,
  t += delt;
  hatu = (hatu + seig*FourierDCT[u^3 - (1 + a)*u])/cheig;
  u = FourierDCT[hatu];
  ], ArrayPlot[Sign[u - Mean[u]] + 1, PixelConstrained -> 1]]
ArrayPlot[Sign[u - Mean[u]] + 1, PixelConstrained -> 1]

After 10000 time steps:

Answer 3

Here's a version using CellularAutomaton to simulate the Ising model (with Glauber dynamics), which can be viewed as a lattice version of the Cahn–Hilliard equation. It's not quite as smooth as the GIF, but it can probably be improved with some tweaks:

With[{L = 100, T = 0.1, t = 100, dt = 10, r = 0.1},
    ArrayPlot[#, ColorRules -> {-1 -> Black, 1 -> White}, PixelConstrained -> True] & /@
    CellularAutomaton[
        {
            ( {{_, a_, _},{d_, x_, b_},{_, c_, _}} ) :> 
                If[RandomReal[] < r/(E^(2/T x (a + b + c + d)) + 1), -x, x]
        },
        RandomChoice[{-1, 1}, {L, L}],
        {{0, t, dt}}]
] // ListAnimate

Answer 4

init = RandomReal[1, {200, 200}];
alpha = -1; beta = -1; gamma = 1; s = 50;

replenish[n_] := RandomReal[1, {n}]
step[a_] := 
  MapThread[
   Prepend, {MapThread[
     Append, {Insert[
       Table[a[[i, j]] + 
         alpha (a[[i, j + 1]] - a[[i, j]] + a[[i, j - 1]]) + 
         beta (a[[i + 1, j + 1]] + a[[i - 1, j + 1]] + 
            a[[i + 1, j - 1]] + a[[i - 1, j - 1]] - a[[i, j]]) + 
         gamma (a[[i + 1, j]] + a[[i - 1, j]] - a[[i, j]]), {i, 2, 
         Length@a - 1}, {j, 2, Length@First@a - 1}], 
       replenish[Length@a - 2], {{-1}, {1}}], replenish[Length@a]}], 
    replenish[Length@a]}];
Do[Subscript[mat, 1] = init; 
 Subscript[mat, i] = step[Subscript[mat, i - 1]];, {i, 2, s}]
ListAnimate[
 Table[ReliefPlot[Subscript[mat, i], ColorFunction -> GrayLevel], {i, 
   s}]]

---

Upon vaxquis's remark in the comments, I am posting a much more concise version of the code. All credit goes to him.

conv = {{γ, α, γ}, {β, 1-α-β-γ, β}, {γ, α, γ}}; 
α = 1; β = -1; γ = -1; 

list = ListConvolve[conv, RandomReal[1, {200, 200}], 2]; 
lists={list}; 
For[i=0, i<50, i++,
  list = ListConvolve[conv, list, 2]; 
  AppendTo[lists, list];
]; 

ListAnimate[ReliefPlot[#, ColorFunction -> GrayLevel]& /@ lists]

---

J.M. suggested the use of NestList[] instead of For[] in the code above:

conv = {{γ, α, γ}, {β, 1-α-β-γ, β}, {γ, α, γ}}; 
α = 1; β = -1; γ = -1; 
ListAnimate[
 ReliefPlot[#, ColorFunction -> GrayLevel] & /@ 
  NestList[ListConvolve[conv, #, 2] &, 
   ListConvolve[conv, RandomReal[1, {200, 200}], 2], 50]]

This decreases the computation time from about 2.96 to 2.82 seconds on my machine, averaged over 10 trials.