Error In Jacobi Iterative Method

I am trying to do part (a) but am having troubles with my code. My solution error is ~300% because my b matrix (the updated u_h(x_i,y_j) matrix) is not calculating properly.

It ends on b =

         0         0         0         0
         0   -1.1111   -1.3333    0.1111
         0   -1.3333   -1.4444    0.4444
         0    0.1111    0.4444    1.0000

where the actual solution matrix is T =

         0         0         0         0
         0    0.0123    0.0494    0.1111
         0    0.0494    0.1975    0.4444
         0    0.1111    0.4444    1.0000

I am not sure where I am making my mistake and I've been trying to fix it for hours now, any help is appreciated. Attached is my code:

 function WIP( maxit,n,eps )
h = 1/(n-1);                % Grid Spacing
x = (0:h:1);                % Define x vector
y = (0:h:1);                % Define y vector
tol = eps;                  % Tolerance
k = 0;                      % Iteration counter
a = zeros(n,n);             % Define residual size
b = zeros(n,n);             % Define solution size
T = zeros(n,n);             % Define true solution size
f = @(X,Y) X^2*Y^2;         % Function
g = @(X,Y) 2*(X^2+Y^2);     % Second Derivative of Function                
done = false;               % Boolean for Iterative Loop
c = 0;                      % Define the constant c

 for k = 1:n                         %Define Boundary Values
        a(1,k) = f(x(1),y(k));
        a(n,k) = f(x(n),y(k));
        a(k,1) = f(x(k),y(1));
        a(k,n) = f(x(k),y(n));
        b(1,k) = f(x(1),y(k));
        b(n,k) = f(x(n),y(k));
        b(k,1) = f(x(k),y(1));
        b(k,n) = f(x(k),y(n));
 end

 for i = 1:n                         %True Solution
    for j = 1:n
        T(i,j) = f(x(i),y(j));
    end
 end

 for i = 1:n                         %INITIAL VALUE OF APPROXIMATION
    for j = 1:n
        b(i,j) = (1-x(i))*a(1,j) + x(i)*a(n,j) + (1-y(j))*a(i,1) + ...
            y(j)*a(i,n) - (1-y(j))*(1-x(i))*a(1,1) - (1-y(j))*x(i)*a(n,1) ...
            - y(j)*(1-x(i))*a(1,n) - x(i)*y(j)*a(n,n);
    end
 end

 while ~done                         %While Loop To Solve Poisson 2D Unit Square
    denom = norm((b-a),inf);        %Difference in solution before Jacobi
    k = k+1;                        %Increase Iteration Counter
    a = b;                          %Update a matrix every iteration

    for j=2:n-1                     %Jacobi Method
        for i=2:n-1
            b(i,j) = (a(i+1,j) + a(i,j+1) + a(i-1,j) + a(i,j-1) - (h^2)*g(i,j))/4;
        end
    end

    numer = norm((b-a),inf);        %Difference in solution after Jacobi

    if c < (numer/denom) && (numer/denom) < 1   %Constant C calculation
       c = numer/denom;
    end

    if numer < tol                  %End While Loop If tolerance is reached
        done = true;
    end;

    if k > maxit                    %End While Loop If Max Iteration is reached
        done = true;
    end

 end

                                    %Print Important Variables
fprintf('\n\nGAUSS-JACOBI METHOD\n\n');
fprintf('N = %d\n\n',n);
fprintf('Number of Iterations = %d\n\n',k);
fprintf('maximum c = %.16f\n\n',c);
est_error = (c/(1-c))*norm((b-a),inf);
true_error = norm((T-b),inf);
fprintf('Estimated Error = %.16f',est_error);
fprintf('\n\nTrue Error = %.16f',true_error);
fprintf('\n\n');

 end

Answer 1

[Torsten]

Your evaluation of the function g is not correct.

Use

g(x(i),y(j))

instead of

g(i,j).

Best wishes

Torsten.