MATLBA官方给出的2D Allen-Cahn Matlab代码分析

PDE

%Solving 2D Allen-Cahn Eq using pseudo-spectral with Implicit/Explicit
%u_t= u_{xx}+u_{yy} + u - u^3 （偏微分方程）
%where u-u^3 is treated explicitly and u_{xx} and u_{yy} is treated implicitly
%BC = Periodic （周期性边界条件）
%IC=v=sin(2*pi*x)+0.001*cos(16*pi*x; （初值条件）

clear all; clc; %清除工作空间的所有变量，函数，和MEX文件。清空命令窗口中的内容。

%Grid （绘制网格）
N = 64; h = 1/N; x = h*(1:N);  %N是网格点个数，h是每个网格大小，x是网格点的空间序列
dt = .01; %dt是时间采样频谱

%x and y meshgrid
y=x'; %y是x的转置
[xx,yy]=meshgrid(x,y); % 表示区域网格设置,目的是为了让x,y成为矩阵格点

%initial conditions （初值条件）
v=sin(2*pi*xx)+0.001*cos(16*pi*xx);
epsilon=.01; %epsilon的值

%(ik) and (ik)^2 vectors in x and y direction
kx=(1i*[0:N/2-1 0 -N/2+1:-1]); %(0+1i)*[0:31 0 -31:-1]
ky=(1i*[0:N/2-1 0 -N/2+1:-1]');%kx的转置
k2x=kx.^2; %kx的平方
k2y=ky.^2; %ky的平方

ii=1:N; %ii=[1 2 3 ...64]
%(ik)^2 matricies in x and y direction
for m= 1:N
    uxx(m,:)= k2x(ii); %Second derivative in the x-direction （x 方向的u的二阶导数）
end
 
for j= 1:N 
    uyy(:,j)= k2y(ii); %Second derivative in the y-direction  （y 方向的u的二阶导数）
end  

for n = 1:500 %时间步长为0.01，跑500次时间步
    
    v_nl=v.^3;  %calculates nonlinear term in real space
    
    for m=1:N %FFT in x-direction on linear and nonlinear term
        v_nl(m,:) = fft(v_nl(m,:)); %对v_n的第m行进行快速Fourier变换
        v_hat(m,:)=fft(v(m,:)); %对v的第m行进行快速Fourier变换
    end
    
    for j=1:N %FFT in y-direction on linear and nonlinear term
        v_nl(:,j) = fft(v_nl(:,j)); %对v_nl的第j列进行快速Fourier变换
        v_hat(:,j)=fft(v_hat(:,j)); %对v_hat的第j列进行快速Fourier变换
    end
    
    vnew=(v_hat/dt-v_nl)./ ...
       (-(uxx+uyy)*epsilon+1/dt-1); %Implicit/Explicit timestepping（隐式/显式时间步长）
  
    for m=1:N %converts to real space in x-direction
    %跑完FFT后我们就得到了多项式乘积的点值表示，现在我们需要将点值表示转回系数表示，这个转换的过程被称为离散傅里叶逆变换（IDFT）
        v(m,:)=ifft(vnew(m,:)); %使用快速傅里叶变换算法计算 vnew的第m行的逆离散傅里叶变换
    end
    
    for j=1:N %converts to real space in y-direction
        v(:,j)=real(ifft(v(:,j)));  %返回对v的第j列进行逆变换的复数的实部数值
    end
    
   %Plots each timestep %画图
   surf(v); 
   title(num2str(n)); 
   axis([0 N 0 N -1 1]); 
   view(43,22); 
   drawnow;
     
end