5
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ I can get the exact solution to get those initial and boundary conditions $\endgroup$ – Mikey Spivak Apr 8 at 14:36
  • $\begingroup$ Just for the sake of seeing the numerical, cant we just use the exact solution to obtain initial and boundary? $\endgroup$ – Mikey Spivak Apr 8 at 18:42
9
+250
$\begingroup$

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;

output

$\endgroup$
  • $\begingroup$ can you show the implementation in matlab of the crank-nicolson? $\endgroup$ – James Apr 29 at 8:21
  • $\begingroup$ Is it possible you can show us the ADI implementation as well? Thanks $\endgroup$ – ILoveMath May 1 at 18:12
  • $\begingroup$ @ILoveMath Yes, I'll do that (needs to modify the matrices) $\endgroup$ – Harry49 May 1 at 18:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.