nx = 41;
ny = 41;
Lx = 1.0;
Ly = 1.0;
A = Lx;
dx = Lx / (nx - 1);
dy = Ly / (ny - 1);
omega = 0.5; 
tol = 1e-5;
max_iter = 5000;
% Initialize Psi, zeta, and temperature T
Psi = zeros(nx, ny); % Initial guess for Psi (stream function)
zeta = zeros(nx, ny); % Initial guess for vorticity
T = zeros(nx, ny); % Temperature
% Boundary conditions for initial time tau = 0
T(:, :) = 0; % Initially, T = 0 throughout
zeta(:, :) = 0; % Initially, zeta = 0 throughout
% Boundary conditions for tau > 0
% X = 0 and X = A (left and right boundaries)
Psi(:, 1) = 0; % Psi = 0 at X = 0
Psi(:, end) = 0; % Psi = 0 at X = A
% dPsi/dX = 0 at X = 0 and X = A (already satisfied as Psi = 0 at boundaries)
% dT/dX = 0 at X = 0 and X = A (no-flux boundary)
% Y = 0 (bottom boundary)
Psi(1, :) = 0; % dPsi/dY = 0 at Y = 0
T(1, :) = -1; % T = -1 at Y = 0
% Y = 1 (top boundary)
Psi(end, :) = 0; % dPsi/dY = 0 at Y = 1
T(end, :) = 1; % T = +1 at Y = 1
% Main SLOR iteration loop
for iter = 1:max_iter
Psi_old = Psi; % Save old Psi for comparison
for i = 2:nx-1
for j = 2:ny-1
% Compute the velocity components U and V using centered differences
U = (Psi(i, j+1) - Psi(i, j-1)) / (2 * dy); % U = dPsi/dY
V = -(Psi(i+1, j) - Psi(i-1, j)) / (2 * dx); % V = -dPsi/dX
% Successive Line Over Relaxation formula for Psi
Psi(i, j) = (1 - omega) * Psi(i, j) + ...
omega / 4 * (Psi(i+1, j) + Psi(i-1, j) + ...
Psi(i, j+1) + Psi(i, j-1) + dx^2 * zeta(i, j));
end
end
% Compute the vorticity zeta using centered differences after Psi update
for i = 2:nx-1
for j = 2:ny-1
% zeta = dV/dX - dU/dY
zeta(i, j) = ((Psi(i+1, j) - Psi(i-1, j)) / (2 * dx)) - ...
((Psi(i, j+1) - Psi(i, j-1)) / (2 * dy));
end
end
% Check for convergence using the relative error
error = max(max(abs(Psi - Psi_old)));
if error < tol
fprintf('Converged after %d iterations with error %e\n', iter, error);
break;
end
end
% If maximum iterations are reached without convergence
if iter == max_iter
disp('Maximum iterations reached without convergence.');
end
% Plotting the stream function Psi
figure;
contourf(linspace(0, Lx, nx), linspace(0, Ly, ny), Psi', 20);
colorbar;
title('Stream Function \Psi using SLOR');
xlabel('x');
ylabel('y');
