solve discrete poisson equation with an input vorticity field using sor in matlab in matlab

main.m
% Define domain parameters
N = 100;        % Number of grid points
L = 1;          % Domain size

% Generate the vorticity field as an example
omega = rand(N);    % Replace with your own vorticity field generation code
219 chars
7 lines

main.m
% Set up grid and spacing
dx = L/N;
x = linspace(0, L, N);
y = linspace(0, L, N);

% Initialize the stream function field
psi = zeros(N);

% Initialize convergence criteria
max_iter = 1000;             % Maximum number of iterations
tolerance = 1e-6;            % Tolerance for convergence
converged = false;           % Convergence flag
iter = 0;                    % Iteration counter
387 chars
14 lines

main.m
% SOR parameters
omega_sor = 1.8;             % Over-relaxation factor

while ~converged && iter < max_iter
    % Perform one SOR iteration
    iter = iter + 1;
    
    for i = 2:N-1
        for j = 2:N-1
            psi(i, j) = (1-omega_sor)*psi(i, j) + omega_sor/(2/dx^2 + 2/dx^2)*(psi(i+1, j) + psi(i-1, j) + psi(i, j+1) + psi(i, j-1) - dx^2*omega(i, j));
        end
    end
    
    % Check convergence
    residual = max(abs(dx^2*(psi(2:N-1, 2:N-1) - (psi(1:N-2, 2:N-1) + psi(3:N, 2:N-1) + psi(2:N-1, 1:N-2) + psi(2:N-1, 3:N-1)) + dx^2*omega(2:N-1, 2:N-1)))));
    if residual < tolerance
        converged = true;
    end
end
634 chars
20 lines

main.m
% Plot the stream function field
figure;
contourf(x, y, psi, 'LineStyle', 'none');
colorbar;
xlabel('x');
ylabel('y');
title('Stream Function Field');

% Add vector plot of velocity field (optional)
[X, Y] = meshgrid(x, y);
u = zeros(N);
v = zeros(N);
for i = 2:N-1
    for j = 2:N-1
        u(i, j) = -(psi(i+1, j) - psi(i-1, j))/(2*dx);
        v(i, j) = (psi(i, j+1) - psi(i, j-1))/(2*dx);
    end
end
hold on;
quiver(X, Y, u, v, 'k');
hold off;
449 chars
22 lines