[Math] Implementing numerical scheme for 2D heat equation in MATLAB

Question

We have 2D heat equation of the form

$$ v_t = \frac{1}{2-x^2-y^2} (v_{xx}+v_{yy}), \; \; \; \; (x,y) \in (-1/2,1/2) \times (-1/2,1/2) $$

We can solve this equation for example using separation of variables and we obtain exact solution

$$ v(x,y,t) = e^{-t} e^{-(x^2+y^2)/2} $$

Im trying to implement the Crank-nicolson and the Peaceman-Rachford ADI scheme for this problem using MATLAB. Can someone help me out how can we do this using matlab?

Answer 1

Let us consider the equation $v_t = \frac{1}{2-x^2-y^2} (v_{xx}+v_{yy})$ with initial and boundary conditions $$v|_{t=0} = e^{-(x^2+y^2)/2}, \quad v|_{y=\pm1/2} = e^{-(t+x^2/2+1/8)}, \quad v|_{x=\pm1/2} = e^{-(t+y^2/2+1/8)}.$$ For the Crank-Nicolson part, we implement the scheme $$ \frac{v_{i,j}^{n+1} - v_{i,j}^n}{k} = \frac{1}{2-{x_i}^2-{y_j}^2} \frac{1}{2 h^2}\left((\delta^2_x + \delta^2_y)v_{i,j}^{n} + (\delta^2_x + \delta^2_y)v_{i,j}^{n+1}\right) , $$ where $v_{i,j}^n \simeq v (x_i,y_j, n k)$ denotes the numerical solution, $(x_i,y_j) = \big(i-\frac{N+1}{2}, j-\frac{N+1}{2}\big) h$ are the grid-node coordinates, $h=\frac{1}{N-1}$ is the mesh size, and $k$ is the time step. The second-order finite difference operators are defined by $$\delta^2_x v_{i,j}^{n} = v_{i+1,j}^n-2v_{i,j}^n+v_{i-1,j}^n, \qquad \delta^2_y v_{i,j}^{n} = v_{i,j+1}^n-2v_{i,j}^n+v_{i,j-1}^n\, .$$ Separating implicit and explicit parts, we rewrite the method as $$ \left(1 - \frac{\lambda_{i,j}}{2} (\delta^2_x + \delta^2_y)\right) v_{i,j}^{n+1} = \left(1 + \frac{\lambda_{i,j}}{2} (\delta^2_x + \delta^2_y)\right) v_{i,j}^{n}, \qquad \lambda_{i,j} = \frac{k/h^2}{2-{x_i}^2-{y_j}^2} $$ for $1< i,j < N$, where $\lambda_{i,j}$ is the mesh Fourier number. Sometimes, the modified method $$ \left(1 - \frac{\lambda_{i,j}}{2} \delta^2_x\right) \left(1 - \frac{\lambda_{i,j}}{2} \delta^2_y\right) v_{i,j}^{n+1} = \left(1 + \frac{\lambda_{i,j}}{2} \delta^2_x\right) \left(1 + \frac{\lambda_{i,j}}{2} \delta^2_y\right) v_{i,j}^{n} $$ is preferred, which approximates the original Crank-Nicolson method up to negligible terms. This method rewrites as the two-step ADI method, which intermediate step ${\bf v}^*$ satisfies \begin{aligned} \left(1 - \frac{\lambda_{i,j}}{2} \delta^2_x\right) v_{i,j}^{*} &= \left(1 + \frac{\lambda_{i,j}}{2} \delta^2_y\right) v_{i,j}^{n} \\ \left(1 - \frac{\lambda_{i,j}}{2} \delta^2_y\right) v_{i,j}^{n+1} &= \left(1 + \frac{\lambda_{i,j}}{2} \delta^2_x\right) v_{i,j}^{*} \, . \end{aligned} A Matlab implementation of Crank-Nicolson is shown below, along with numerical results.

% mesh settings
nx = 20;
x = linspace(-0.5,0.5,nx);
y = linspace(-0.5,0.5,nx);
Fo = 0.5; % k/h^2
tmax = 5e-2;
h = max(abs(diff(x)));
k = Fo*h^2;

% initialization
[Xi,Yi] = meshgrid(x(2:end-1),y(2:end-1));
X = reshape(Xi,(nx-2)^2,1);
Y = reshape(Yi,(nx-2)^2,1);
Uinit = exp(-0.5*(Xi.^2+Yi.^2));
U = Uinit;
u = reshape(U,(nx-2)^2,1);
Lambdai = Fo./(2-Xi.^2-Yi.^2);
lambdai = reshape(Lambdai,(nx-2)^2,1);
t = 0;

% graphics
figure(1);
subplot(1,2,1);
psurf = surf(Xi,Yi,U);
ptit = title(strcat('t =',num2str(t)));
xlim([-0.5,0.5]);
ylim([-0.5,0.5]);
zlim([0,1]);
subplot(1,2,2);
pnum = plot(Xi(floor((nx-1)/2),:),U(floor((nx-1)/2),:),'b.'); hold on
pthe = plot(Xi(floor((nx-1)/2),:),exp(-t)*Uinit(floor((nx-1)/2),:),'k-');
xlim([-0.5,0.5]);
ylim([0,1]);

% matrices
dm = 0.5*lambdai(2:end).*(abs(Y(2:end)+0.5-h)>eps);
dp = 0.5*lambdai(1:end-1).*(abs(Y(1:end-1)-0.5+h)>eps);
A = diag(1+2*lambdai) - diag(0.5*lambdai(nx-1:end),-(nx-2)) - diag(0.5*lambdai(1:end-nx+2),nx-2) - diag(dm,-1) - diag(dp,1);
A = sparse(A);
b = -0.5*lambdai.*(exp(-0.5*X.^2).*((abs(Y+0.5-h)<=eps)+(abs(Y-0.5+h)<=eps)) + exp(-0.5*Y.^2).*((abs(X+0.5-h)<=eps)+(abs(X-0.5+h)<=eps)))*exp(-0.125);
M = 2*inv(A) - eye((nx-2)^2);

% Crank-Nicolson iterations
while (t+k<tmax)
    u = M*(u - b*exp(-t)) - b*exp(-(t+k));
    t = t+k;

    U = reshape(u,nx-2,nx-2);
    set(ptit,'String',strcat('t =',num2str(t)));
    set(pnum,'YData',U(floor((nx-1)/2),:));
    set(pthe,'YData',exp(-t)*Uinit(floor((nx-1)/2),:));
    drawnow;
end

set(psurf,'ZData',U);
drawnow;