# Multigrid method
By 李恩齊 & 姚博宸
[toc]
## Matlab Code
For a evenly-spaced $N\times N$ matrix $\textbf{U}$ that satisfy the following Poisson's equation:
\begin{align}
\frac{\partial u(x,y)}{\partial x}+\frac{\partial u(x,y)}{\partial y}=f(x,y),
\end{align} where the grid is squared, i.e., $\Delta x=\Delta y$, $\textbf{U}$ can be solved with the following lines of code:
### The main function

function U_next = multigrid(U,F,dx)
    s = 3;
    for r = 1:s
        U = gaussseidel(U,F,dx);
    end
    R  = residual(U,F,dx);
    Rc = R(1:2:end,1:2:end);
    if length(Rc)-1 > 2
        Ec = multigrid(zeros(size(Rc)),Rc,2*dx);
    else        
        Ec = gaussseidel(zeros(size(Rc)),Rc,2*dx);
    end
    E = interp_c2f(U,Ec);
    U = U-E;
    for r = 1:s
        U = gaussseidel(U,F,dx);
    end
    U_next = U;
end

### Gauss-Seidel method

function U_next = gaussseidel(U,F,dx)
    U_next = U;
    for i = 2:size(U,1)-1
        for j = 2:size(U,1)-1
            U_next(i,j) = (U(i-1,j)+U(i+1,j)+U_next(i,j-1)+U_next(i,j+1)-dx^2*F(i,j))/4;
        end
    end
end

### The residual matrix

function R = residual(U,F,dx)
    R = zeros(size(U));
    R(2:end-1,2:end-1) = (U(1:end-2,2:end-1)+U(3:end,2:end-1)+U(2:end-1,1:end-2)+U(2:end-1,3:end)-4*U(2:end-1,2:end-1))/dx^2-F(2:end-1,2:end-1);
end

### Interpolation from coarse grid to fine grid

function E = interp_c2f(U,Ec)
    E       = zeros(size(U));
    [Xc,Yc] = meshgrid(linspace(0,1,length(Ec)));
    [X,Y]   = meshgrid(linspace(0,1,length(E)));
    E       = interp2(Xc,Yc,Ec,X,Y);
end