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

Question

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} = 2Delta 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),; tin(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)

I love math · Answer

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’)

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

One Answer

Add your own answers!

Ask a Question