MATLAB code: Third-order Runge-Kutta Fourier spectral method for phase-field crystal model with FCC ordering structure

%2D code for phase transition

Lx = 256; xleft = 0; xright = xleft+Lx;

Ly = 256; yleft = 0; yright = yleft+Ly;

Mx = 256; dx = Lx/Mx; x = xleft+(0:Mx-1)*dx;

My = 256; dy = Ly/My; y = yleft+(0:My-1)*dy;

xix = 2pi[0:Mx/2-1 -Mx/2:-1]/Lx;

xiy = 2pi[0:My/2-1 -My/2:-1]/Ly;

[kx,ky] = ndgrid(xix,xiy);

dt = 0.1;

h = dx;

R = 0.01(2rand(Mx,My)-1); ophi = 0.23+(R-mean(mean(R)));

epsilon = 0.195; g = 0; r = 0.0; s = 2; T = 2000; order = 3;

qq = sqrt(2); rr = 0.01;

kp = -(s-(2qq^4+2qq^2+2rr^2)(kx.^2+ky.^2)+…

(qq^4+4qq^2+1+rr^2)(kx.^2+ky.^2).^2-…

(2qq^2+2)(kx.^2+ky.^2).^3+(kx.^2+ky.^2).^4);

switch order case 1

MI = 1; ME = MI;

case 2

gamma = (2-sqrt(2))/2; lap = 1-1/(2*gamma);

MI = [gamma 0; 1-gamma gamma]; ME = [gamma 0; lap 1-lap];

case 3

MI = [1/2 0 0 0; 1/6 1/2 0 0; -1/2 1/2 1/2 0; 3/2 -3/2 1/2 1/2];

ME = [1/2 0 0 0; 11/18 1/18 0 0; 5/6 -5/6 1/2 0; 1/4 7/4 3/4 -7/4];

end

ns = length(MI);

for n=1:round(T/dt)

ophi_hat = fft2(ophi);

phis = zeros(Mx,My,ns); q = zeros(Mx,My,ns);

for i = 1:ns

Phip = ophi.^3-epsilonophi+(qq^4+rr^2)ophi;

Phir = (1/4ophi.^4-epsilon/2ophi.^2+(qq^4+rr^2)/2*ophi.^2).^r;

beta = sum(sum(Phip))/sum(sum(Phir));

q(:,:,i) = -(Phip-sophi)+Phirbeta;

IM = 0; EX = ME(i,i)*q(:,:,i);

for j=1:i-1

IM = IM + MI(i,j)*phis(:,:,j);

EX = EX + ME(i,j)*q(:,:,j);

end

phis(:,:,i) = (ophi_hat+dt(kp.IM+fft2(EX)))./(1-MI(i,i)dtkp);

ophi = real(ifft2(phis(:,:,i)));

end

if (mod(n,40)==0 )

figure(1);clf;

[xx,yy] = ndgrid(x,y); pcolor(xx,yy,ophi); shading interp;

colormap jet; axis image off;

end

end