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 = Lb1;
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?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);

%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 = Lb1;
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? 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);

%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 = Lb1;
u = u;
u = reshape(u,N,N);

Cart

Cart