PDE solver I made for the laplace equation fails when I insert a poisson equation

Hello all, I wrote a solver to the poisson equation, it actually works fine for the laplace equation (e.g, when the laplacian equals to 0), but it screws by a lot when the laplacian isn't equal to 0.

function Sol = LaplaceSolver(xLimit, yLimit, boundaryFunction, rectFunction, N_x, N_y,tol, LinearSolverType)
    %% in this function we solve the Laplace Equation using Jacobi's method
    %% first we define our grid
    
    x_s = linspace(xLimit(1), xLimit(2), N_x);

    y_s = linspace(yLimit(1), yLimit(2), N_y);

    delta_x = x_s(2)-x_s(1);
    delta_y = y_s(2)-y_s(1);

    %% We essentially solve the linear system of equation Ax = b where:
    %% b = boundary
    %% A the matrix
    %% we begin with the bottom layer

    b = DefineBoundaryLaplace5(x_s,y_s, boundaryFunction, rectFunction);

    %% Next we determine our system of equations
    delta_x_sq = delta_x^2;
    delta_y_sq = delta_y^2;
    FiniteDifferenceEquationLaplace = DefineEquationLaplace4(N_x,N_y,delta_x_sq,delta_y_sq);

    %% Finally we use the Jacobi's method to solve our (sparse) linear system of equation
    xInitial = zeros(length(b),1);
    switch LinearSolverType
        case "Jacobi"
        VectorSol = JacobiMethod(FiniteDifferenceEquationLaplace,b,xInitial, tol);
        case "Gauss"
        VectorSol = GaussSiedelMethod(FiniteDifferenceEquationLaplace,b,xInitial, tol);
        case "SOR"
        omega = 4 / (2 + sqrt(4 - (cos(pi / N_y) + cos(pi/N_x)^2)));
        VectorSol = SORMethod(FiniteDifferenceEquationLaplace,b,xInitial,omega,tol);
        case "ConjGrad"
            VectorSol = ConjGrad(FiniteDifferenceEquationLaplace,b,tol);
        case "weird"
            VectorSol = Weirdy(b,tol)
    end
    %% and trully finally, we reshape our VectorSol into a matrix:
    Sol = reshape(VectorSol, [N_x-2,N_y-2]);


end

The function DefineBoundaryLaplace5:

function b = DefineBoundaryLaplace5(x,y, boundaryFunction, rectFunction)
    N_x = length(x);
    N_y = length(y);
    
    delta_x = x(2) - x(1);
    delta_y = x(2) - x(1);
    b = zeros((N_x-2),(N_y-2));
    
    delta_x_sq = delta_x^2;
    delta_y_sq = delta_y^2;

    arr_x = 1:(N_x-2);
    arr_y = 1:(N_y-2);

    x_coeffs = 1;
    y_coeffs = delta_x_sq / delta_y_sq;
    rect_coeffs = delta_x_sq;

    b(arr_x,arr_y) = rect_coeffs * rectFunction(x(arr_x+1), y(arr_y+1));
    b(arr_x,1) = b(arr_x,1) + y_coeffs * boundaryFunction(x(arr_x+1), y(1))';
    b(arr_x, N_y-2) = b(arr_x,N_y-2) + y_coeffs * boundaryFunction(x(arr_x+1), y(N_y))';
    b(1, arr_y) = b(1,arr_y) + x_coeffs * boundaryFunction(x(1), y(arr_y+1));
    b(N_x-2, arr_y) = b(N_x-2,arr_y) + x_coeffs * boundaryFunction(x(N_x), y(arr_y+1));

    b = reshape(b, [(N_x-2)*(N_y-2),1]);
end

The function DefineEquationLaplace4:

function FiniteDifferenceLaplace = DefineEquationLaplace4(N_x,N_y,delta_x_sq,delta_y_sq)
    FiniteDifferenceLaplace = cell((N_x-2) * (N_y-2),1);

    L = CalcL(1,1,N_y);

    
    x_coeff = -1;
    y_coeff = -delta_x_sq / delta_y_sq;

    Center_coeff = 2 * (delta_x_sq / delta_y_sq + 1);
    FiniteDifferenceLaplace{L} = [Center_coeff,1, x_coeff,2, y_coeff,N_x-1]; %% bottom left
    L = CalcL(N_x-2,1,N_y);
    FiniteDifferenceLaplace{L} = [x_coeff,L-1,Center_coeff,L,y_coeff,L+(N_x-2)]; % bottom right
    L = CalcL(1, N_y-2, N_y);
    FiniteDifferenceLaplace{L} = [y_coeff,L-(N_x-2),Center_coeff,L,x_coeff,L+1]; % upper left
    L = CalcL(N_x-2, N_y-2, N_y);
    FiniteDifferenceLaplace{L} = [y_coeff,L-(N_x-2),x_coeff,L-1,Center_coeff,L]; % upper right


    for i = 2:(N_x-3) % i iterates on x
        L = CalcL(i,1,N_y); % bottom line
        FiniteDifferenceLaplace{L} = [x_coeff,L-1,Center_coeff,L,x_coeff,L+1,y_coeff, L + (N_x-2)]; % bottom line
        L = CalcL(i,N_y-2,N_y); % upper line
         FiniteDifferenceLaplace{L} = [y_coeff, L - (N_x-2),x_coeff,L-1,Center_coeff,L,x_coeff,L+1]; % upper line

    end

     for j = 2:(N_x-3) % j iterates on y 
        L = CalcL(1,j,N_y); % left column
        FiniteDifferenceLaplace{L} = [y_coeff,L-(N_x-2),Center_coeff,L,x_coeff,L+1,y_coeff, L + (N_x-2)]; % left line

        L = CalcL(N_x-2,j,N_y); % upper line
        FiniteDifferenceLaplace{L} = [y_coeff,L-(N_x-2),x_coeff,L-1,Center_coeff,L,y_coeff, L + (N_x-2)]; % right line
    end


    for i = 2:(N_x-3) % i iterates on x
        for j = 2:(N_y-3) % j iterates on y
            L = CalcL(i,j,N_y);
            FiniteDifferenceLaplace{L} = [y_coeff,L - (N_x-2),x_coeff,L-1,Center_coeff,L,x_coeff,L+1,y_coeff, L + (N_x-2)]; % all other circumstances
        end
    end
    
end

end

And for instance, the Jacobi method which I use to iteratively solve the equation Ax = b:

function x = JacobiMethod(Matrix,b, xInitial, tol)
%% this function solves a sparse linear system of equations
%% it receives a cell array A{1},...,A{n}, each cell represents a line in
%% the matrix in the following way - A{1} = {A(1i1, i1),i1, A(1i2),i2, A1(i3),i3),...}
%% this is an efficent way to represent a sparse matrix, since we simply dont represent the zeros.
%% A = [[a1_i1, i1 , a_i2, i2, ], [ai_i, i,...]


N = length(xInitial);
x = zeros(N,1);
for i = 1:N
    CurrLine = Matrix{i};
    
    SumDotProduct = sum(CurrLine(1:2:length(CurrLine)) * xInitial(CurrLine(2:2:length(CurrLine))) );


    diagonal_id = find(CurrLine(2:2:length(CurrLine)) == i);
    

    a_ii = CurrLine(diagonal_id * 2 - 1);
    SumDotProduct = SumDotProduct - a_ii * xInitial(i);
    x(i) = ( -SumDotProduct + b(i)) / a_ii;
end


norm_to_check = max(abs(x - xInitial)) / max(abs(x));

if (norm_to_check < tol)
    return
else
x = JacobiMethod(Matrix,b,x,tol);
end

end

Thanks in advance for anyone willing to help, and I will also appericate any other advises on how to make my code even more optimized (even though it runs pretty smooth when I use the Conj grad)

Answer 1

[UDAYA PEDDIRAJU]

Hi Matan,

To address the issue with poison equation solver you can try the following:

1. Ensure Proper Source Term Integration: Check that the source term is correctly added to the “b” vector.
2. Check Boundary Conditions: Confirm that “DefineBoundaryLaplace5” applies the boundary conditions accurately for both equations.
3. Verify Matrix Coefficients: Make sure that “DefineEquationLaplace4” accurately represents the discretized Poisson equation.
4. Review Convergence Criteria: Adjust the tolerance and iteration count for the Poisson equation as needed.
5. Test Solver Stability: Experiment with different initial guesses and consider using a more robust solver if necessary.
6. Optimize Code: Use vectorization and efficient sparse matrix operations to enhance performance.
7. Debugging: Use debugging tools to examine intermediate computations, especially after incorporating the source term.
8. Assess Numerical Stability: Check if the discretization scheme requires modification for stability with the source term.

Additionally, you can refer to the following toolbox, that can solve the Poisson equation: https://www.mathworks.com/matlabcentral/fileexchange/38090-2d-poisson-equation.

I hope this helps!