How to write a code of 2D ADI method in matlab?

I tried to write a code for the alternating direction implicit (ADI) method in 2D, but I got stuck.

My equation is:

$$\frac{\partial U(t,x,y)}{\partial t} = 2\Delta U(t,x,y) -10(\frac{\partial U(t,x,y)}{\partial x} + \frac{\partial U(t,x,y)}{\partial x}) + [(x^2 + y^2 + 4t(5x + 5y -2)] \; \forall \;(x,y) \in (0,1)\times(0,1),\; t\in(0,T]$$

Initial and boundary conditions are:

\begin{align} &U(0,x,y) = 0 \;\forall\;(x,y)\in[0,1]\times[0,1]\\ &U(t,x,y) =t (x^2 + y^2) \; \forall \; (x,y) \in \partial[0,1]\times[0,1] \end{align}

I solved this equation (I used ADI method) and I found:

1.step $$(-15)u_{i-1,j}^{k+1/2} + (26)u_{i,j}^{k+1/2} + (-10)u_{i+1,j}^{k+1/2} = (15)u_{i,-1j}^{k} + (-24)u_{i,j}^{k} + (10)u_{i,j+1}^{k} + [\frac{x_i^2}{2} + \frac{y_j^2}{2} + 2t(5x + 5y -2)]$$ from Dirichlet condition: $$i=1 \qquad \qquad u_{0,j}^{k+1/2}= t_{k+1/2}(0^2 + y_j^2) \\ i=n-1 \; \qquad u_{n,j}^{k+1/2}=t_{k+1/2}(1^2+y_j^2)$$

2.step $$(-15)u_{i-1,j}^{k+1} + (26)u_{i,j}^{k+1} + (-10)u_{i+1,j}^{k+1} = (15)u_{i,-1j}^{k+1/2} + (-24)u_{i,j}^{k+1/2} + (10)u_{i,j+1}^{k+1/2} + [\frac{x_i^2}{2} + \frac{y_j^2}{2} + 2t(5x + 5y -2)]$$ from Dirichlet condition: $$j=1 \qquad \qquad u_{i,0}^{k+1}= t_{k+1}(x_i^2 + 0^2) \\ j=n-1 \; \qquad u_{i,n}^{k+1}=t_{k+1}(x_i^2+1^2)$$

Although I found both steps, but I do not know, hot to write it in Matlab. Could somebody help me, please?

My attempt is:

n = 10;
h = 1/n;
tau = 1/10;
mi = tau/h/h;

g0x=@(t,y)0.*x;
g1x=@(t,y)0.*x;
g0y=@(t,y)0.*y;
g1y=@(t,y)0.*y;

u0=@(x,y)0,*+0,*y;

x = 0 : h : 1;
y = 0 : h : 1;
t = 0 : tau:1;
[X,Y]=meshgrid(x,y)

u = zeros(J+1,Nt);

for i = 1:n+1
u(i,n)=(-15)*u(i-1,n) + 26*u(i,n) + (-10)*u(i+1,n)+ (-15)*u(i,n-1) + 24*u(i,n)+(-10)u(i,n+1)
end
u(0,n)=t*u(y*y)
u(n,n)=t*u(1+y*y)

for j = 1:n+1
u(n,j)=(-15)*u(n,j-1) + 26*u(n,j) + (-10)*u(n,j+1)+ (-15)*u(n,j-1) + 24*u(n,j)+(-10)u(n,j+1)
end
u(0,n)=t*u(y*y)
u(n,n)=t*u(1+y*y)

• @nicoguaro I think that I understand, how to use ADI, but I do not know, how to write it in Matlab. Could you help me with this code, please? You do not "have to" write a whole code ... there is possibility that I just need some hint ... I wrote equations for the 1st and 2nd step with conditions, but I do not know, how to write it in Matlab ... – I love math May 25 at 14:59

I found a code for ADI in 1D, but I do not know, how to change it and I do not understand why I need an analytical solution of my PDE (this code is in book Computational Partial Differential Equations Using MATLAB - Jichun Li, Yi-Tung Chen). Mabye it could help somebody to understand, what I need.

% para1d.m:
% para1d.m:
% use the explicit scheme to solve the parabolic equation
% u_t(x,t) = u_{xx}(x,t), xl < x < xr, 0 < t < tf
% u(x,0) = f(x), xl < x < xr
% u(0,t) = gl(t), u(1,t) = gr(t), 0 < t < tf
%
% A special case is choosing f and g properly such that the
% analytic solution is:
% u(x,t)= sin(pi*x)*e^(-pi^2*t) + sin(2*pi*x)*e^(-4*pi^2*t)
% we solve this program by the explicit scheme:
% u(j,n+1) = u(j,n) + v*(u(j+1,n) - 2*u(j,n) + u(j-1,n))
%---------------------------------------------------------------
clear all; % clear all variables in memory
xl=0; xr=1; % x domain [xl,xr]
J = 10; % J: number of division for x
dx = (xr-xl) / J; % dx: mesh size
tf = 0.1; % final simulation time
Nt = 50; % Nt: number of time steps
dt = tf/Nt;
mu = dt/(dx)^2;
if mu > 0.5 % make sure dt satisy stability condition
error(’mu should < 0.5!’)
end
% Evaluate the initial conditions
x = xl : dx : xr; % generate the grid point
% f(1:J+1) since array index starts from 1
f = sin(pi*x) + sin(2*pi*x);
% store the solution at all grid points for all time steps
u = zeros(J+1,Nt);
% Find the approximate solution at each time step
for n = 1:Nt
t = n*dt; % current time
% boundary condition at left side
gl = sin(pi*xl)*exp(-pi*pi*t)+sin(2*pi*xl)*exp(-4*pi*pi*t);
% boundary condition at right side
gr = sin(pi*xr)*exp(-pi*pi*t)+sin(2*pi*xr)*exp(-4*pi*pi*t);
if n==1 % first time step
for j=2:J % interior nodes
u(j,n) = f(j) + mu*(f(j+1)-2*f(j)+f(j-1));
end
u(1,n) = gl; % the left-end point
u(J+1,n) = gr; % the right-end point
else
for j=2:J % interior nodes
u(j,n)=u(j,n-1)+mu*(u(j+1,n-1)-2*u(j,n-1)+u(j-1,n-1));
end
u(1,n) = gl; % the left-end point
u(J+1,n) = gr; % the right-end point
end
% calculate the analytic solution
for j=1:J+1
xj = xl + (j-1)*dx;
u_ex(j,n)=sin(pi*xj)*exp(-pi*pi*t) ...
+sin(2*pi*xj)*exp(-4*pi*pi*t);
end
end
% Plot the results
tt = dt : dt : Nt*dt;
figure(1)
colormap(gray); % draw gray figure
surf(x,tt, u’); % 3-D surface plot
xlabel(’x’)
ylabel(’t’)
zlabel(’u’)
title(’Numerical solution of 1-D parabolic equation’)
figure(2)
surf(x,tt, u_ex’); % 3-D surface plot
xlabel(’x’)
ylabel(’t’)
zlabel(’u’)
title(’Analytic solution of 1-D parabolic equation’)

• This does not seem as an answer to the OP. – nicoguaro May 23 at 14:06
• I found it in book which I mentioned above and I thought that it could help somebody to answer my question .... – I love math May 25 at 15:01