1
$\begingroup$

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)
$\endgroup$
  • $\begingroup$ @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 ... $\endgroup$ – I love math May 25 at 14:59
0
$\begingroup$

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’)
| cite | improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ This does not seem as an answer to the OP. $\endgroup$ – nicoguaro May 23 at 14:06
  • $\begingroup$ I found it in book which I mentioned above and I thought that it could help somebody to answer my question .... $\endgroup$ – I love math May 25 at 15:01

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.