# Modified diffusion equation and unstabilities

I am trying to simulate the phase separation of a binary mixture. If the free energy F is known as a function of the concentration $c$, the dynamical equation is:

$\frac{\partial c(x,t)}{\partial t}=\frac{d^2}{dx^2} \frac{\delta F[c]}{\delta c}$

For the Flory-Huggins free energy we have:

$\frac{\delta F[c]}{\delta c}=\log\left(\frac{c}{1-c}\right) + \chi c^2 + \gamma (c')^2$

First term is entropy, second term is attraction between particles ($\chi<0$) and third term is similar as a surface tension.

I use time forward and space centered differentiation. Even for very small $dt$ I get numerical instabilities. I first thought the term $\gamma(c'^2)$ was responsible, but instabilities remain even without.

Here is my Matlab code for no-flux boundary conditions, clear and simplified as much as I can. I am aware the boundary condition implementation may not be correct but I don't think the problem comes from this.

What should I do?

function phaseSep ()
clear all;
clf;

N = 101;   %number of grid points in x
dx = 1/(N - 1);
x = 0:dx:1;   %vector of x values
T = 1e3;  %number of time steps
dt = 1e-8;

% Second derivative
function y=mylaplace(f,i)
y = f(i+1)-2*f(i)+f(i-1);
end

function y=dfdc(C)

p=-0.01;
g=0.01;

for i=2:N-1
y(i) = log(C(i)/(1-C(i)))+p*C(i)+g*mylaplace(C,i)/dx^2;
end

% No-flux Boundary conditions
y(1)=y(2);
y(N)=y(N-1);
end

%Initial concentration
for i=1:N
C(i)=0.2+0.1*tanh(10*(x(i)-0.5));
end

% Plot initial concentration
hold;
plot(x,C);

% iterate
for t=1:T
Z=dfdc(C);
for i=2:N-1
C(i) = C(i) + (mylaplace(Z,i)/dx^2)*dt;
end

%No flux boundary conditions
C(1) = C(2);
C(N) = C(N-1);
end

% Plot final concentration
plot(x,C);

end

• What happens if you take the simplest possible flux, $F[c]=c^2/2$? – Wolfgang Bangerth Jan 26 '15 at 2:52
• So this reduce to the diffusion equation and that works fine for $dt \leq 10^{-5}$ – David Jan 26 '15 at 12:56
• @geoff put my thoughts in words. An alternative venue to get an explanation of why you have to choose time steps so small is in lectures 26 and 27 at math.tamu.edu/~bangerth/videos.html – Wolfgang Bangerth Jan 27 '15 at 4:06

At first glance, it looks like you are using the method of lines with forward Euler time steps. What WolfgangBangerth is getting at with his example is that even for a simple heat equation, the stability limit of forward Euler (namely, that $|\lambda\Delta{t}| < 1$) combined with the eigenvalues induced by a finite difference approximation of the diffusion operator result in a very stringent limit on the time step relative to the spatial discretization. For a second-order centered finite difference approximation to the second-order derivative term, this limit is of the form $\Delta{t} < h^{2}/(2a)$, where $a$ is the diffusivity and $h$ is the grid spacing of the spatial discretization (assumed uniform). You should check out a text such as LeVeque's Finite Difference Methods for Ordinary and Partial Differential Equations for a basic overview of the theory. In terms of practical advice, you probably want to replace the time discretization you are currently using with an L-stable implicit method. To start, you might try using ode23tb in MATLAB.