Jacobi and Gauss-Seidel Method for PDE

Figure 1 shows both analytical and numerical solution contours for a 41x41 mesh. The numerical solution was performed by Jacobi method.

Figure 2 shows the residual (R2) as a function of iteration for 41x41 and 81x81 mesh for both Jacobi and Gauss-Seidel method.

%%(Jacobi Method)

% Arif Hossain

clc; close all; clearvars;

mesh = 81;

N = mesh;

x = linspace(0,1,N);

y = linspace(0,1,N);

[xx,yy] = meshgrid(x,y);

M = N;

L = 1;

H = 1;

delx = L/(N-1);

dely = H/(M-1);

%% Analytical calculation

for j = 1:M

    for i = 1:N

        phi_anly(j,i) = 500*exp(-50*((1-x(i)).^2+y(j).^2))+100*x(i)*(1-y(j));

    end

end

%% Numerical calculation

% Boundary condition

phi(1:M,1:N) = 0; % Initial guess

tolerance = 1e-6;

A = 1/(delx.^2);

B = 1/(dely.^2);

C = 2*(A+B);

count = 0; 

R22 = 0;

for j = 1:M

    for i = 1:N

        phi(1,i)= 100*x(i)+500*exp(-50*(1-x(i))^2);

        phi(N,i)= 500*exp(-50*((1-x(i))^2+1));

        phi(j,1)= 500*exp(-50*(1+y(j)^2));

        phi(j,N)= 100*(1-y(j))+500*exp(-50*y(j)^2);

        p(j,i) = exp(-50*((1-x(i))^2+y(j)^2));

        q(j,i) = 100*((1-x(i))^2+y(j)^2)-2;

        s(j,i) = 50000*p(j,i)*q(j,i);

    end

end

%% Jacobi Method

for iter = 1:100000

    count = count+1;

    phi_star = phi;

    for j = 2:M-1

        for i = 2:N-1

phi(j,i) = ((A*phi_star(j,i+1)+A*phi_star(j,i-1)+B*phi_star(j 1,i)+B*phi_star(j+1,i))- s(j,i))/C;

        end

    end 

    clear i j

    R2 = 0;

    for j = 2:M-1

        for i = 2:N-1

            R(j,i) = s(j,i)-A*phi(j,i+1)-A*phi(j,i-1)-B*phi(j-1,i)-B*phi(j+1,i)+C*phi(j,i);

            R2 = R2+R(j,i)*R(j,i);

        end

    end

    R2 = sqrt(R2);

    if R2<tolerance

        break

    end  

R22(count) = R2;

end   

%% Residual plot

    figure()

    semilogy(R22,'-r','linewidth',3)

    set(gca,'fontsize',14,'fontweight','bold')

    xlabel('Iteration'),ylabel('Residual,R2');grid on

    ylim([1e-6 1e6])

%% Error   

for j = 2:N

    for i = 2:M  

        phi_err(j,i) = (((phi_anly(j,i)-phi(j,i))/(phi_anly(j,i))))*100;

    end

end

%% Contour Plot

        figure()

        contourf(xx,yy,phi,15)

        set(gca,'fontsize',14,'fontweight','bold')

        xlabel('x'),ylabel('y');grid on

        title('\phi_n (Numerical)');

        colormap(jet);colorbar

        figure()

        contourf(xx(1:N-1,2:M),yy(1:N-1,2:M),phi_err(1:N-1,2:M),15)

        set(gca,'fontsize',14,'fontweight','bold')

        xlabel('x'),ylabel('y');grid on

        title('Error (%)');

        colormap(jet);colorbar

        caxis([0 1])

%%(Gauss-Seidel Method)

% Arif Hossain

clc; close all; clear all;

mesh = 81;

N = mesh;

M = N;

L = 1;

H = 1;

delx = L/(N-1);

dely = H/(M-1);

x = linspace(0,1,N);

y = linspace(0,1,N);

[xx,yy] = meshgrid(x,y);

%% Analytical calculation

for j = 1:M

    for i = 1:N

        phi_anly(i,j) = 500*exp(-50*((1-x(j)).^2+y(i).^2))+100*x(j)*(1-y(i));

    end

end

%% Numerical calculation

% Boundary condition

phi(1:M,1:N) = 0; % Initial guess

tolerance = 1e-6;

A = 1/(delx.^2);

B = 1/(dely.^2);

C = 2*(A+B);

count = 0; 

R22 = 0;

for j = 1:M

    for i = 1:N

        phi(1,i)= 100*x(i)+500*exp(-50*(1-x(i))^2);

        phi(N,i)= 500*exp(-50*((1-x(i))^2+1));

        phi(j,1)= 500*exp(-50*(1+y(j)^2));

        phi(j,N)= 100*(1-y(j))+500*exp(-50*y(j)^2);

        p(j,i) = exp(-50*((1-x(i))^2+y(j)^2));

        q(j,i) = 100*((1-x(i))^2+y(j)^2)-2;

        s(j,i) = 50000*p(j,i)*q(j,i); % source term

    end

end

%% Gause-Seidel Method

for iter = 1:100000

    count = count+1;

    for j = 2:M-1

        for i = 2:N-1       

              phi(j,i) = (A*phi(j,i+1)+A*phi(j,i-1)+B*phi(j-1,i)+B*phi(j+1,i)-s(j,i))/C;

        end

    end

    R2 = 0;

    for j = 2:M-1

        for i = 2:N-1

            R(j,i) = s(j,i)-A*phi(j,i+1)-A*phi(j,i-1)-B*phi(j-1,i)-B*phi(j+1,i)+C*phi(j,i);

            R2 = R2+R(j,i)*R(j,i);

        end

    end

    R2 = sqrt(R2);

    if R2<tolerance

        break

    end  

R22(count) = R2;

end   

%% Residual plot

    figure()

    semilogy(R22,'-r','linewidth',3)

    set(gca,'fontsize',14,'fontweight','bold')

    xlabel('Iteration'),ylabel('Residual,R2');grid on

    ylim([1e-6 1e6])

%% Error   

for j = 2:N

    for i = 2:M  

        phi_err(j,i) = (((phi_anly(j,i)-phi(j,i))/(phi_anly(j,i))))*100;

    end

end

%% Contour Plot

        figure()

        contourf(xx,yy,phi,15)

        set(gca,'fontsize',14,'fontweight','bold')

        xlabel('x'),ylabel('y');grid on

        title('\phi_n (Numerical)');

        colormap(jet);colorbar

        figure()

        contourf(xx(1:N-1,2:M),yy(1:N-1,2:M),phi_err(1:N-1,2:M),15)

        set(gca,'fontsize',14,'fontweight','bold')

        xlabel('x'),ylabel('y');grid on

        title('Error (%)');

        colormap(jet);colorbar

        caxis([0 1])