Incorporate the SOR Method into the multigrid test case MATLAB code.

SOR MATLAB CODE:

% SOR (Successive Over-Relaxation)

n = input('Enter number of equations, n: ');

A = zeros(n,n+1);

x1 = zeros(1,n);

A=[5 -2 3 -1;-3 9 1 2;2 -1 -7 3];

x1 = [0 0 0];

tol = input('Enter tolerance, tol: ');

m = input('Enter maximum number of iterations, m: ');

w = input('Enter the parameter w (omega): ');

k = 1;

while k <= m

err = 0;

for i = 1 : n

s = 0;

for j = 1 : n

s = s-A(i,j)*x1(j);

end

s = w*(s+A(i,n+1))/A(i,i);

if abs(s) > err

err = abs(s);

end

x1(i) = x1(i)+s;

end

if err <= tol

break;

else

k = k+1;

end

end

fprintf('The solution vector after %d iterations is : ', k);

for i = 1 : n

fprintf(' %11.8f ', x1(i));

end

MULTIGRID TEST CASE MATLAB CODE:

%% Multigrid test case with visulization

% Author: Ben Beaudry

clc; clear; close all

load A_b.mat;

A = full(A); % Unsparce matrix to show full power of Multigrid

pre = 2; % Number of presmoothing iterations

post = 2; % Number of postsmoothing iterations

% cycle = 1; % Type of multigrid cycle (1=V-cycle, 2=W-cycle, 3=F-cycle)

smooth = 1; % Smoother type (1=Jacobi, 2=Gauss-Seidel)

grids = 5; % Max grids in cycle

maxit = 10; % Max iterations of solver

tol = 1e-02; % Tolerance of solver

%% Solvers

% solve directly

t=cputime;

x_D = A\b;

fprintf('Direct Solve Time: %g Seconds ',cputime-t)

% V-Cycle

t=cputime;

[x_V,res_V,it_V] = multigrid(A,b,pre,post,1,smooth,grids,maxit,tol);

fprintf('V-Cycle Solve Time: %g Seconds ',cputime-t)

% W-Cycle

t=cputime;

[x_W,res_W,it_W] = multigrid(A,b,pre,post,2,smooth,grids,maxit,tol);

fprintf('W-Cycle Solve Time: %g Seconds ',cputime-t)

% F-Cycle

t=cputime;

[x_F,res_F,it_F] = multigrid(A,b,pre,post,3,smooth,grids,maxit,tol);

fprintf('F-Cycle Solve Time: %g Seconds ',cputime-t)

% max it for iterative methods

it = 100;

% Gauss-Siedel

t=cputime;

L = tril(A);

U = triu(A,1);

x_G = zeros(length(b),1);

res_G= zeros(it,1);

for g = 1:it

x_G = L\(b-U*x_G);

r = b - A * x_G;

res_G(g) = norm(r);

if res_G(g) < tol

break;

end

end

fprintf('Gauss-Seidel Solve Time: %g Seconds ',cputime-t)

% Jacobi

t=cputime;

d = diag(A);

dinv = 1./d;

LU = tril(A,-1)+triu(A,1);

U = triu(A,1);

x_J = zeros(length(b),1);

res_J= zeros(it,1);

for j = 1:it

x_J = dinv.*(b-(LU*x_J));

r = b - A * x_J;

res_J(j) = norm(r);

if res_J(j) < tol

break;

end

end

fprintf('Jacobi Solve Time: %g Seconds ',cputime-t)

%% Plotting

figure;

hold on

plot(1:g,res_G(1:g),'o-','DisplayName','Guass-Seidel')

plot(1:j,res_J(1:j),'o-','DisplayName','Jacobi')

plot(1:it_V,res_V(1:it_V),'o-','DisplayName','V-Cycle Multigrid')

plot(1:it_F,res_F(1:it_F),'o-','DisplayName','F-Cycle Multigrid')

plot(1:it_W,res_W(1:it_W),'o-','DisplayName','W-Cycle Multigrid')

set(gca,'XLim',[0 100]);

grid on

legend('location','best')

ylabel('Relative Convergance')

xlabel('Iterations')

title('Convergance of Numerical System')

hold off