Clear Filters
Clear Filters

Need help solving heat equation using adi method

11 views (last 30 days)
Vard
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
% ADI method
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;
  1 Comment
Torsten
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.

Sign in to comment.

Answers (1)

John D'Errico
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...
  19 Comments
Vard
Vard on 13 May 2024
@Torsten Do you have a another solution with using three dimensional arrays?
Torsten
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

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!