Need help solving heat equation using adi method
6 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
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 Commento
Torsten
il 12 Mag 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.
Risposte (1)
John D'Errico
il 12 Mag 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...
Vedere anche
Categorie
Scopri di più su Eigenvalue Problems in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!