Modified diffusion equation and unstabilities

Question

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

Answer 1

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.