MATLAB code: Fourier spectral method for Molecular Beam Epitaxy (MBE) model

clear all;

Nx= 128; Ny= 128;

Lx=2pi ; Ly=2pi ;

Lx=64 ; Ly=64 ;

hx=Lx/Nx ; hy=Ly/Ny;

x=linspace(-0.5Lx+hx,0.5Lx,Nx);

y=linspace(-0.5Ly+hy,0.5Ly,Ny) ;

epsilon = 0.5;

[xx,yy]= ndgrid(x,y);

%u = 0.1( sin ( 3 *xx ) .sin ( 2 * yy ) +sin ( 5 xx ) .sin ( 5 * yy ) ) ;

for i = 1:Nx

for j = 1:Ny

     u(i,j) = 0.001*(2.0*rand-1);

end

end

p=1i*2 * pi/Lx *[ 0 :Nx/2-1 0 -Nx/2+1:-1];

q=1i*2 * pi/Ly *[ 0 :Ny/2-1 0 -Ny/2+1:-1] ;

[pp,qq]= ndgrid(p,q);

p2 = ( 2* pi/Lx *[ 0 :Nx/2 -Nx/ 2 + 1 : - 1 ] ).^ 2 ;

q2 = ( 2* pi/Ly *[ 0 :Ny/2 -Ny/ 2 + 1 : - 1 ] ).^ 2 ;

[pp2,qq2]= ndgrid(p2,q2);

dt = 0.001; T=145; Nt=round ( T/dt ) ; ns=Nt/50;

figure(1) ; clf; colormap gray;

surf(x,y,real(u’)); axis image ; view(0,90); shading interp;

%colorbar ;

caxis ( [ -1 1 ] ) ; pause(0.01);

for iter = 1:Nt

u=real(u); tu= fft2(u);

fx=real(ifft2(pp.*tu));

fy=real(ifft2(qq.*tu));

f1 = (fx.^2 + fy.^2).*fx ;

f2 = (fx.^2 + fy.^2).*fy ;

s_hat = fft2(u/dt) + (pp.fft2(f1) + qq .fft2(f2)) ;

v_hat = s_hat./( 1/dt - ( pp2+qq2 ) + epsilon*( pp2+qq2 ).^2 ) ;

u= ifft2( v_hat ) ;

if (mod( iter , ns ) ==0)

figure(2); clf ;

subplot(1,2,1); colormap gray; axis image ;

surf( x, y, real(u’)) ; view ( 0,90 ) ; shading interp ;

subplot(1,2,2); colormap gray; axis image ;

lapu = ifft2(-(pp2+qq2).*v_hat);

surf( x, y, real(lapu’)) ; view ( 0,90 ) ; shading interp ;

caxis ( [ -1 1 ] ) ;

pause (0.01)

end

end