create a multigrid solver for poissons equation in matlab

main.m
% Create a multigrid solver for Poisson's equation in MATLAB

% Define the parameters
L = 1; % Length of the domain
N = 64; % Number of grid points
h = L / (N + 1); % Grid spacing
max_levels = log2(N) - 1; % Maximum number of multigrid levels

% Create the finest grid
x = linspace(0, L, N+2); % Grid points
f = -1 * ones(1, N+2); % Right-hand side
u = zeros(1, N+2); % Initial guess

% Define the relaxation parameter and number of iterations
omega = 2/3; % Over-relaxation parameter
num_iters = 10; % Number of relaxation iterations

% Create multigrid levels
levels = cell(max_levels, 1);
levels{1} = struct('u', u, 'f', f, 'h', h, 'x', x);

for k = 2:max_levels
    Nk = 2^(k-1) * N;
    hk = L / (Nk + 1);
    xk = linspace(0, L, Nk+2);
    fk = -1 * ones(1, Nk+2);
    levels{k} = struct('u', zeros(1, Nk+2), 'f', fk, 'h', hk, 'x', xk);
end

% Multigrid cycle
for cycle = 1:10
    % V-cycle
    for k = 1:max_levels
        u = vcycle(levels, u, k, omega, num_iters);
    end
end

% Function for the V-cycle
function u = vcycle(levels, u, k, omega, num_iters)
    if k == 1
        u = relax(u, levels{k}.f, omega, num_iters);
    else
        u = relax(u, levels{k}.f, omega, num_iters);
        r = residual(u, levels{k}.f, levels{k}.h);
        r = restrict(r);
        e = zeros(size(r));
        e = vcycle(levels, e, k-1, omega, num_iters); % Recursion
        e = interpolate(e);
        u = u + e;
        u = relax(u, levels{k}.f, omega, num_iters);
    end
end

% Relaxation
function u = relax(u, f, omega, num_iters)
    for iter = 1:num_iters
        for i = 2:length(u)-1
            u(i) = (1-omega)*u(i) + (omega/2) * (u(i-1) + u(i+1) + h^2*f(i));
        end
    end
end

% Residual calculation
function r = residual(u, f, h)
    r = zeros(size(u));
    for i = 2:length(u)-1
        r(i) = f(i) - (u(i-1) - 2*u(i) + u(i+1)) / h^2;
    end
end

% Restriction
function rc = restrict(r)
    rc = zeros(1, length(r)/2 + 1);
    for i = 1:length(rc)
        rc(i) = r(2*i-1) + 0.25*(r(2*i) + r(2*i-2));
    end
end

% Interpolation
function ef = interpolate(e)
    ef = zeros(1, 2*length(e)-1);
    for i = 1:length(e)
        ef(2*i) = e(i);
        ef(2*i-1) = 0.5 * (e(i) + e(i+1));
    end
end
2216 chars
87 lines