MATLAB code: Energy-stable linear convex splitting method for the Quasi-crystal pattern

clear all;

Lx = 128; xleft = 0; xright = xleft+Lx;

Ly = 128; 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 = 2^(-3);h = dx;R = 0.1(2rand(Mx,My)-1); ophi = 0.0+(R-mean(mean(R)));

epsilon = 0.015;alpha = 1.0; T = 512;qq = 2*cos(pi/12); beta = 1.0;

gamma = -epsilon + 3beta(2*alpha/3+beta);

kp = -(qq^4-(2qq^4+2qq^2)(kx.^2+ky.^2)+(qq^4+4qq^2+1)*(kx.^2+ky.^2).^2…

-(2qq^2+2)(kx.^2+ky.^2).^3+(kx.^2+ky.^2).^4);Ip = zeros(Mx,My);phis = zeros(Mx,My);

for n=1:round(T/dt)

ophi_hat = fft2(ophi);

%first stage:

for i = 1:Mx

for j = 1:My

if(ophi(i,j) > 2*alpha/3 + beta)

Ip(i,j) = 3beta(2alpha/3 + beta)ophi(i,j) - (alpha/3+2beta)(2*alpha/3+beta)^2;

elseif(ophi(i,j) >= -beta && ophi(i,j) <= 2*alpha/3+beta)

Ip(i,j) = ophi(i,j)^3 - alpha*ophi(i,j)^2;

else

Ip(i,j) = 3beta(2alpha/3+beta)ophi(i,j) + beta^2(alpha+2beta);

end

end

end

source1 = ophi_hat - dt( -(gamma+epsilon)ophi_hat + fft2(Ip) );

phis = source1./(1+dtgamma-dtkp);ophi1 = real(ifft2(phis));

%second stage:

ophi_hat1 = fft2(ophi1);

for i = 1:Mx

for j = 1:My

if(ophi1(i,j) > 2*alpha/3 + beta)

Ip(i,j) = 3beta(2alpha/3 + beta)ophi1(i,j) - (alpha/3+2beta)(2*alpha/3+beta)^2;

elseif(ophi1(i,j) >= -beta && ophi1(i,j) <= 2*alpha/3+beta)

Ip(i,j) = ophi1(i,j)^3 - alpha*ophi1(i,j)^2;

else

Ip(i,j) = 3beta(2alpha/3+beta)ophi1(i,j) + beta^2(alpha+2beta);

end

end

end

source2 = -0.5ophi_hat + 1.5ophi_hat1 - 0.5dt( -(gamma+epsilon)*ophi_hat1 + fft2(Ip) );

phis = source2./(1+0.5dtgamma-0.5dtkp);ophi2 = real(ifft2(phis));

%third stage:

ophi_hat2 = fft2(ophi2);

for i = 1:Mx

for j = 1:My

if(ophi2(i,j) > 2*alpha/3 + beta)

Ip(i,j) = 3beta(2alpha/3 + beta)ophi2(i,j) - (alpha/3+2beta)(2*alpha/3+beta)^2;

elseif(ophi2(i,j) >= -beta && ophi2(i,j) <= 2*alpha/3+beta)

Ip(i,j) = ophi2(i,j)^3 - alpha*ophi2(i,j)^2;

else

Ip(i,j) = 3beta(2alpha/3+beta)ophi2(i,j) + beta^2(alpha+2beta);

end

end

end

source3 = -0.5ophi_hat + 2.5ophi_hat1 - ophi_hat2…

• 0.5dt( -(gamma+epsilon)*ophi_hat2 + fft2(Ip) );

phis = source3./(1+0.5dtgamma-0.5dtkp);ophi = real(ifft2(phis));

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

figure(1);clf;

[xx,yy] = ndgrid(x,y); pcolor(xx,yy,ophi); shading interp; colormap jet; axis image off;

end

end