Kronecker product implementation of finite difference for poisson equation

Hi all,

I want to test the setup to use kron to solve the 2D poisson equation with dirichlet boundary condition on two sides (left and bottom) and non-zero neumann boundary conditions on the other sides. I know the theoretical solution which is and thus the forcing function is . I have the code herein:

%grid declaration 
N = 100;
x = linspace(0,1,N);
dx = x(2)-x(1);
y = x;
[X,Y] = meshgrid(x,y);

%analytical solution and forcing function 
u_an = 8*X.^2.*Y.^2;
f = 16*(X.^2+Y.^2);

%creation of the matrix (complete matrix including end points)
e = ones(N,1);
T = spdiags([e -2*e e],[-1 0 1],N,N);
T(1,1:2) = [dx^2 0]; %dirichlet matrix (first order approximation)
T(end,end-1:end) = [dx -dx]; %neumann matrix (first order approximation)
fullT = full(T);
I = speye(N);
L = kron(T,I)+kron(I,T);
L = L/dx^2;
fullL = full(L);

%rhs including boundary conditions and forcing function 
b = zeros(N,N);
b(2:end-1,2:end-1) = f(2:end-1,2:end-1);
b(:,end) = -16*y'.^2; %right boundary
b(end,:) = -16*x.^2; %top boundary

%solution of the matrix equation and reshaping 
b1 = reshape(b,N^2,1);
u = L\b1;
u = u; 
u = reshape(u,N,N);

%visualisation and comparison with the analytical solution 
figure(1)
plot3(X,Y,u,'*k')
hold on 
plot3(X,y,u_an,'or')
hold off

The code isn't quite accurate but what is bizzare is that if you replace

T(end,end-1:end) = [dx -dx]; %neumann matrix (first order approximation)

and

b(:,end) = -16*y'.^2; %right boundary
b(end,:) = -16*x.^2; %top boundary

with

T(end,end-1:end) = [1 -1]; %neumann matrix (first order approximation)

and

b(:,end) = -16*y'.^2/dx; %right boundary
b(end,:) = -16*x.^2/dx; %top boundary

it generates the correct result. Could someone please help me understand if I'm setting up the Kronecker product correctly or not and what changes in the two implementations?

[Torsten]

Maybe of interest:

https://uk.mathworks.com/matlabcentral/fileexchange/27279-laplacian-in-1d-2d-or-3d?s_tid=FX_rc1_behav

Answer 1

[Umar]

Hi @Aditya,

@Torsten provided file exchange link which is not a bad idea to check. You probably know this that in order to solve the 2D Poisson equation using Kronecker products effectively, it's essential to ensure that both the matrix representation of the differential operator and the right-hand side (RHS) vector correctly reflect the specified boundary conditions. Let me help you understand matrix setup which will help clarify confusion.

Understanding Matrix Setup

1. Matrix (T): The matrix (T) represents a one-dimensional discrete Laplacian operator, which approximates the second derivative. For Dirichlet conditions, you set (T(1,1:2) = [dx^2, 0] at the left edge, which correctly implements a fixed value (Dirichlet). Now, for Neumann conditions, your original setup ( T(end,end-1:end) = [dx, -dx] ) tries to enforce a zero-gradient (Neumann), but it does so incorrectly by using (dx) which scales with grid spacing rather than enforcing a constant gradient.

2. Boundary Condition Vector ( b ): The right-hand side vector (b) must reflect your boundary conditions accurately. When you implemented (b(:,end) = -16*y'.^2; )and (b(end,:) = -16*x.^2; ), these terms were not scaled correctly for Neumann conditions since they need to reflect flux rather than values at boundaries. You can double check this.

Here is my modified approach:

T(end,end-1:end) = [1 -1]; % Correct for Neumann condition
b(:,end) = -16*y'.^2/dx; % Correct scaling for right boundary
b(end,:) = -16*x.^2/dx; % Correct scaling for top boundary

This change will ensure that:

*The Neumann condition is enforced correctly by setting it up as a finite difference approximation for derivatives (using differences rather than scaled values).

*The RHS incorporates (dx) scaling to reflect that these are derivatives (fluxes), which is appropriate for Neumann boundaries.

Here are some additional points to help further

Kronecker Product Usage: Your use of `L = kron(T,I) + kron(I,T); is appropriate as it constructs a 2D Laplacian from a 1D version effectively. Just ensure that both components correctly reflect your boundary conditions.

Numerical Stability: Ensure that your grid size (N)is sufficiently large to minimize numerical errors due to discretization.

Visualization: After solving, comparing your numerical solution against an analytical solution (as you did with u_an) is an excellent practice for validating your implementation.

Further Testing: It might be beneficial to test various grid sizes and compare results to confirm stability across different configurations.

In nutshell, your initial confusion stemmed from misapplying Neumann boundary conditions in both matrix construction and RHS definition. By correcting these elements, you aligned your numerical approach with theoretical expectations, yielding accurate results.

Hope this helps.

[Aditya]

Thanks a lot for clarification! Really helpful!