Effacer les filtres
Effacer les filtres

Need help solving heat equation using adi method

18 vues (au cours des 30 derniers jours)
Vard
Vard le 12 Mai 2024
Modifié(e) : Vard le 16 Mai 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 commentaire
Torsten
Torsten le 12 Mai 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.

Connectez-vous pour commenter.

Réponses (1)

John D'Errico
John D'Errico le 12 Mai 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 commentaires
Vard
Vard le 13 Mai 2024
@Torsten Do you have a another solution with using three dimensional arrays?
Torsten
Torsten le 13 Mai 2024
Modifié(e) : Torsten le 13 Mai 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

Connectez-vous pour commenter.

Community Treasure Hunt

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

Start Hunting!

Translated by