# How to make a Cahn–Hilliard model GIF

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["http://i.stack.imgur.com/y9np9.gif"]

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


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

So how to make this GIF by its first frame?

-
So those images came from this site: en.m.wikipedia.org/wiki/Cahn%E2%80%93Hilliard_equation so you need to translate those equations describing a phase transition into a code – JasonB Feb 14 at 7:16
@JasonB Oh,It looks more difficult than what I think. – yode Feb 14 at 7:34
Please make the question self-contained and improve the title. For instance, Export["/tmp/foo.gif", Import["http://i.stack.imgur.com/y9np9.gif"]] seems to satisfy the requirements, accurately reproducing the gif, and it is quite the easiest way. – Michael E2 Feb 14 at 16:10
@MichaelE2 Sorry my poor English,I have changed some place as your advice just now.If there are other flaw you think in my expression,help to improve it.I will be so thankful. – yode Feb 15 at 17:36
@yode I thought you wanted an implementation of the Cahn-Hilliard model, and I thought you should add a brief description of what it is as a guide to answerers and other readers of the Q&A. It seems (now, to me) that you're after anything that looks similar to the animation, in which case the question seems ok now. – Michael E2 Feb 16 at 13:28

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];


-
Very neat that you can capture the overall effect with image transformations, I never would have thought to try this route. +1 – JasonB Feb 14 at 13:43
Concern the concision and ingeniousness of your second method.I firm accepted your answer in spite of it may not be accurately in accordance with the theorem of Cahn–Hilliard equation .And thanks for your help. – yode Feb 15 at 16:35
@yode Thank you! I'm glad you liked it: It was fun to play with. – MarcoB Feb 15 at 16:44

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:

-
This is the best answer so far IMO, but the code gives a lot of errors and doesn't return a result on my version 10.3 on Windows 10. – thedude Feb 14 at 11:19
@thedude it seems I accidentaly forgot to copy two lines of code so eps2 was not defined. Can you try now? – shrx Feb 14 at 11:24
Now it works fine, quite fast too. – thedude Feb 14 at 11:25
@shrx +1, it's always cool to see interesting physics simulated with Mathematica! – JasonB Feb 14 at 13:43

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


-
Hey Stephen, I took the liberty of adding an animation generated from your code. It's a very cool use of a cellular automaton! (+1) – MarcoB Feb 14 at 21:57
Thanks, @MarcoB ! – Stephen Powell Feb 14 at 22:17
I really like this, well done! – blochwave Feb 15 at 7:51
init = RandomReal[1, {200, 200}];
alpha = -1; beta = -1; gamma = 1; s = 50;

replenish[n_] := RandomReal[1, {n}]
step[a_] :=
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.

-
Where are from the second gif? – yode Feb 14 at 8:06
alpha=-1; beta=1; gamma=-1 – thedude Feb 14 at 8:14
The improvement is so obvious that I wanna give second upvote. – yode Feb 15 at 2:45
Could've used NestList[] instead of a loop: ListAnimate[ReliefPlot[#, ColorFunction -> GrayLevel] & /@ NestList[ListConvolve[conv, #, 2] &, ListConvolve[conv, RandomReal[1, {200, 200}], 2], 50]]. – J. M. Feb 15 at 15:00
@J.M. I can't try it right now, does it make a noticeable performance change? – thedude Feb 15 at 15:50