# Need help solving heat equation using adi method

11 views (last 30 days)
Vard on 12 May 2024
Edited: Vard on 16 May 2024
Nx = 50;
Ny = 50;
dx = 1.0 / (Nx - 1);
dy = 2.0 / (Ny - 1);
T = 0.00001;
dt = 0.000000001;
Nt = 1000; % Total number of time steps
alpha = 4; % Diffusion coefficient
x = linspace(0, 1, Nx);
y = linspace(0, 2, Ny);
[X, Y] = meshgrid(x, y);
U = Y .* X + 1; % Initial condition
function val = f(x, y, t)
val = exp(t) * cos(pi * x / 2) * sin(pi * y / 4);
end
function U = apply_boundary_conditions(U, x, y, dy, dx)
U(:, 1) = 1; % y=0
U(:, end) = U(:, end-1) + dy .* x'; % y=2
U(1, :) = U(2, :) - dx * y; % x=0
U(end, :) = y' + 1; % x=1
end
% Thomas algorithm
function x = thomas_algorithm(a, b, c, d)
n = length(d);
c_star = zeros(n-1, 1);
d_star = zeros(n, 1);
x = zeros(n, 1);
c_star(1) = c(1) / b(1);
d_star(1) = d(1) / b(1);
for i = 2:n-1
temp = b(i) - a(i-1) * c_star(i-1);
c_star(i) = c(i) / temp;
d_star(i) = (d(i) - a(i-1) * d_star(i-1)) / temp;
end
d_star(n) = (d(n) - a(n-1) * d_star(n-1)) / b(n);
x(n) = d_star(n);
for i = n-1:-1:1
x(i) = d_star(i) - c_star(i) * x(i+1);
end
end
for n = 1:Nt
U = apply_boundary_conditions(U, x, y, dy, dx);
% First half-step: X-direction implicit, Y-direction explicit
for j = 2:Ny-1
a = -alpha * dt / (2 * dx^2) * ones(Nx-1, 1);
b = (1 + alpha * dt / dx^2) * ones(Nx, 1);
c = -alpha * dt / (2 * dx^2) * ones(Nx-1, 1);
d = U(j, 2:end-1)' + 0.5 * alpha * dt / dy^2 * (U(j+1, 2:end-1) - 2 * U(j, 2:end-1) + U(j-1, 2:end-1))' + dt * f(x(2:end-1), y(j), n*dt);
U(j, 2:end-1) = thomas_algorithm(a, b(2:end-1), c, d);
end
% Second half-step: Y-direction implicit, X-direction explicit
for i = 2:Nx-1
a = -alpha * dt / (2 * dy^2) * ones(Ny-1, 1);
b = (1 + alpha * dt / dy^2) * ones(Ny, 1);
c = -alpha * dt / (2 * dy^2) * ones(Ny-1, 1);
d = U(2:end-1, i) + 0.5 * alpha * dt / dx^2 * (U(2:end-1, i+1) - 2 * U(2:end-1, i) + U(2:end-1, i-1)) + dt * f(x(i), y(2:end-1), n*dt);
U(2:end-1, i) = thomas_algorithm(a, b(2:end-1), c, d);
end
U = apply_boundary_conditions(U, x, y, dy, dx);
end
% Visualization
surf(X, Y, U);
title('ADI Solution at T=' + string(T));
xlabel('X');
ylabel('Y');
zlabel('U');
colorbar;
Torsten on 12 May 2024
Shouldn't your solution array be three-dimensional instead of two-dimensional ? The way you arranged the code, you only get one solution at time T - the complete history for U is overwritten.

John D'Errico on 12 May 2024
You say it works for sufficiently small values dt.
With that exp(t) term in there, do you seriously expect it to work well for large values of dt? I have dreams myself somedays...
Vard on 13 May 2024
@Torsten Do you have a another solution with using three dimensional arrays?
Torsten on 13 May 2024
Edited: Torsten on 13 May 2024
If your code were correct, you could simply save the solution matrix U after each time step in a three-dimensional matrix:
U_3d = zeros(Nt,Ny,Nx)
for nt = 1:Nt
...
U_3d(nt,:,:) = U;
end