Skip to content

Commit 889ddba

Browse files
mdumettmdumett
committed
addGral combines pure periodic and addBC
addBC allowed periodic BC imposing the conditions explicitly. Per files solve periodic BC without imposing any conditions. The domain does not include boundaries. Gral combines traditional operators and handling periodic boundary conditions without specifying explicitly them
1 parent 0cac0fc commit 889ddba

30 files changed

Lines changed: 1525 additions & 2 deletions

‎examples/matlab/elliptic1DGralNonPerBC.m‎

Lines changed: 34 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,34 @@
1+
% ====================== Test 1 =====================
2+
% 1D Poisson BVP: Dirichlet, Dirichlet Homogeneous BC
3+
% - u'' = 1, 0 < x < 1, u(0) = 0, u(1) = 0
4+
% exact solution: u(x) = x(1-x)/2
5+
% ===================================================
6+
% example that uses lapGral1D
7+
%
8+
close all; clc;
9+
10+
addpath('../../src/matlab');
11+
12+
k = 2;
13+
bvp = 1;
14+
m = 2*k+1;
15+
dx = 1/m;
16+
% centers and vertices
17+
xc = [0 dx/2:dx:1-dx/2 1]';
18+
t = '- u" = 1, 0 < x < 1, u(0) = 0, u(1) = 0, with exact solution u(x) = x(1-x)/2';
19+
ue = 0.5*xc.*(1-xc); % exact solution
20+
dc = [1;1];
21+
nc = [0;0];
22+
v = [0;0];
23+
A = - lapGral1D(k,m,dx,dc,nc);
24+
b = ones(size(A,2),1);
25+
[A0,b0] = addGralBC1D(A,b,k,m,dx,dc,nc,v);
26+
ua = A0\b0; % approximate solution
27+
28+
% plot
29+
figure(bvp)
30+
plot(xc,ue,'b*',xc,ua,'ro');
31+
title(t); %,'interpreter','latex');
32+
xlabel('x');
33+
ylabel('u');
34+
legend({'exact','approx'});

‎examples/matlab/elliptic1DGralPerBC.m‎

Lines changed: 36 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,36 @@
1+
% ====================== Test 5 =====================
2+
% 1D Poisson BVP: Periodic BC
3+
% - u'' = 4 pi^2 sin(2 pi x), 0 < x < 1, u(0) = u(1), u'(0) = u'(1)
4+
% exact solution: u(x) = sin(2 pi x) + constant
5+
% ===================================================
6+
% example that uses lapGral1D
7+
% testing 1D bc
8+
%
9+
close all; clc;
10+
11+
addpath('../../src/matlab');
12+
13+
k = 2;
14+
bvp = 5;
15+
m = 20;
16+
dx = 1/m;
17+
% centers and vertices
18+
xc = (dx/2:dx:1-dx/2)';
19+
t = '- u" = 4 pi^2 sin(2 pi x), 0 < x < 1, u(0) = u(1), u''(0) = u''(1), with exact solution u(x) = sin(2 pi x) + constant';
20+
ue = sin(2*pi*xc); % exact solution
21+
dc = [0;0];
22+
nc = [0;0];
23+
v = [0;0];
24+
A = - lapGral1D(k,m,dx,dc,nc);
25+
b = 4*pi^2 * sin(2*pi*xc);
26+
[A0,b0] = addGralBC1D(A,b,k,m,dx,dc,nc,v);
27+
ua = A0\b0; % approximate solution (there are infinity solutions)
28+
ua = ua - ua(1) + ue(1); % shifting ua to match ue(1) with ua(1)
29+
30+
% plot
31+
figure(bvp)
32+
plot(xc,ue,'b*',xc,ua,'ro');
33+
title(t); %,'interpreter','latex');
34+
xlabel('x');
35+
ylabel('u');
36+
legend({'exact','approx'});

‎examples/matlab/elliptic1DPerBC.m‎

Lines changed: 31 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,31 @@
1+
% ====================== Test 5 =====================
2+
% 1D Poisson BVP: Periodic BC
3+
% - u'' = 4 pi^2 sin(2 pi x), 0 < x < 1, u(0) = u(1), u'(0) = u'(1)
4+
% exact solution: u(x) = sin(2 pi x) + constant
5+
% ===================================================
6+
% example that does not use addBC1D
7+
%
8+
close all; clc;
9+
10+
addpath('../../src/matlab');
11+
12+
k = 2;
13+
bvp = 5;
14+
m = 20;
15+
dx = 1/m;
16+
% centers and vertices
17+
xc = (dx/2:dx:1-dx/2)';
18+
t = '- u" = 4 pi^2 sin(2 pi x), 0 < x < 1, u(0) = u(1), u''(0) = u''(1), with exact solution u(x) = sin(2 pi x) + constant';
19+
ue = sin(2*pi*xc); % exact solution
20+
A = - lapPer(k,m,dx);
21+
b = 4*pi^2 * sin(2*pi*xc);
22+
ua = A\b; % approximate solution (there are infinity solutions)
23+
ua = ua - ua(1) + ue(1); % shifting ua to match ue(1) with ua(1)
24+
25+
% plot
26+
figure(bvp)
27+
plot(xc,ue,'b*',xc,ua,'ro');
28+
title(t); %,'interpreter','latex');
29+
xlabel('x');
30+
ylabel('u');
31+
legend({'exact','approx'});

‎examples/matlab/elliptic2DPer.m‎

Lines changed: 38 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,38 @@
1+
% ====================== Test 2 =====================
2+
% 2D Poisson BVP: Periodic, Periodic domain
3+
% u_xx + u_yy = exp(- 10(x^2 + y^2)), -1 < x,y < 1,
4+
% BC: periodic
5+
% exact solution: unknown
6+
% ===================================================
7+
% example that does not use addBC2D
8+
%
9+
close all; clc;
10+
11+
addpath('../../src/matlab');
12+
13+
k = 2;
14+
bvp = 2;
15+
m = 49; % it should be odd
16+
n = m+2; % it should be odd
17+
dx = 2/m;
18+
dy = 2/n;
19+
% centers and vertices
20+
xc = (-1+dx/2:dx:1-dx/2)';
21+
yc = (-1+dy/2:dy:1-dy/2)';
22+
[Y,X] = meshgrid(yc,xc);
23+
% t = 'u_xx + u_yy = exp(-10(x^2+y^2)), -1 < x,y < 1, periodic boundary conditions. Unknown exact solution';
24+
A = - lap2DPer(k,m,dx,n,dy);
25+
b = - exp(-10*(X.^2 + Y.^2));
26+
b = reshape(b,[],1);
27+
ua = A\b; % approximate solution
28+
ua = reshape(ua,m,n);
29+
ua = ua - ua((m+1)/2,(n+1)/2);
30+
31+
figure(bvp)
32+
surf(X,Y,ua);
33+
title('Approximate Solution: 2D Poisson with Periodic BC');
34+
shading interp;
35+
figure(bvp+10)
36+
surf(X,Y,b);
37+
title('Source term: 2D Poisson with Periodic BC');
38+
shading interp;

‎examples/matlab/elliptic2DXPerYDirichlet.m‎

Lines changed: 46 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,46 @@
1+
% ====================== Test 5 =====================
2+
% 2D Poisson BVP: Periodic BC along X-axis and Dirichlet along Y-axis
3+
% -(u_xx + u_yy) = 2 sin(2 pi x) (1+2 pi^2 y(1-y)), 0 < x,y < 1, u(x,0) = 0 = u(x,1)
4+
% exact solution: u(x) = y(1-y)sin(2 pi x)
5+
% ===================================================
6+
% example that uses addGralBC2D
7+
%
8+
close all; clc;
9+
10+
addpath('../../src/matlab');
11+
12+
k = 2;
13+
bvp = 5;
14+
m = 20;
15+
n = m+1;
16+
dx = 1/m;
17+
dy = 1/n;
18+
% centers and vertices
19+
xc = (dx/2:dx:1-dx/2)';
20+
yc = [0 dy/2:dy:1-dy/2 1]';
21+
[Y,X] = meshgrid(yc,xc);
22+
t = '-(u_xx + u_yy) = 2 sin(2 pi x) (1+2 pi^2 y(1-y)), 0 < x,y < 1, periodic along x, u(x,0) = 0 = u(x,1), with exact solution u(x) = y(1-y)sin(2 pi x)';
23+
dc = [0;0;1;1];
24+
nc = [0;0;0;0];
25+
bcl = 0; % zeros(n,1);
26+
bcr = 0; % zeros(n,1);
27+
bct = zeros(m,1);
28+
bcb = zeros(m,1);
29+
v = {bcl;bcr;bcb;bct};
30+
ue = Y.*(1-Y).*sin(2*pi*X); % exact solution
31+
A = - lapGral2D(k,m,dx,n,dy,dc,nc);
32+
b = 2*sin(2*pi*X).*(1+2*pi^2*Y.*(1-Y));
33+
b = reshape(b,[],1);
34+
[A0, b0] = addGralBC2D(A, b, k, m, dx, n, dy, dc, nc, v);
35+
ua = A0\b0; % approximate solution (there are infinity solutions)
36+
ua = reshape(ua,m,n+2);
37+
38+
% plot
39+
figure(bvp)
40+
surf(X,Y,ua);
41+
title('Approximate Solution: 2D Poisson with Periodic BC along X and Dirichlet on Y');
42+
shading interp;
43+
figure(bvp+10)
44+
surf(X,Y,ue);
45+
title('Exact Solution: 2D Poisson with Periodic BC along X and Dirichlet on Y');
46+
shading interp;

‎examples/matlab/elliptic3DXDirichletYPerZNeumann.m‎

Lines changed: 86 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,86 @@
1+
% ====================== Test 2 =====================
2+
% 3D Poisson BVP: Dirichlet (on left and right faces), Periodic along Y-axis, Neumann along Z-axis BC
3+
% -(u_xx + u_yy + u_zz) = 2 sin(2 pi y) z (1 + 2 pi^2 x(1-x)), 0 < x,y,z < 1
4+
% BC: u(-1,y,z) = 0, u(1,x,y) = 0, periodic along Y-axis, u_z(x,y,0) = -x(1-x)sin(2 pi y), u_z(x,y,0) = x(1-x)sin(2 pi y)
5+
% exact solution: x(1-x)sin(2 pi y)z
6+
% ===================================================
7+
% example that uses addGral3D
8+
%
9+
close all; clc;
10+
11+
addpath('../../src/matlab');
12+
13+
k = 2;
14+
bvp = 2;
15+
m = 29; % should be odd to be able to plot the middle slice
16+
n = m+2; % should be odd to be able to plot the middle slice
17+
o = m+4; % should be odd to be able to plot the middle slice
18+
dx = 1/m;
19+
dy = 1/n;
20+
dz = 1/o;
21+
% centers and vertices
22+
xc = [0 dx/2:dx:1-dx/2 1]';
23+
yc = (dy/2:dy:1-dy/2)';
24+
zc = [0 dz/2:dz:1-dz/2 1]';
25+
[Y,X,Z] = meshgrid(yc,xc,zc);
26+
t = '-(u_xx + u_yy + u_zz) = 2 sin(2 pi y) z (1 + 2 pi^2 x(1-x)), 0 < x,y,z < 1, u(-1,y,z) = 0, u(1,x,y) = 0, periodic along Y-axis, u_z(x,y,0) = -x(1-x)sin(2 pi y), u_z(x,y,0) = x(1-x)sin(2 pi y), with exact solution u(x,y,z) = x(1-x)sin(2 pi y)z';
27+
ue = X.*(1-X).*sin(2*pi*Y).*Z;
28+
dc = [1;1;0;0;0;0];
29+
nc = [0;0;0;0;1;1];
30+
bcl = zeros(n*o,1);
31+
bcr = zeros(n*o,1);
32+
bcb = 0; % zeros((m+2)*o,1);
33+
bct = 0; % zeros((m+2)*o,1);
34+
bcf = zeros(n*(m+2),1);
35+
bcz = zeros(n*(m+2),1);
36+
% bcf = zeros((n+2)*(m+2),1);
37+
% bcz = zeros((n+2)*(m+2),1);
38+
v = {bcl;bcr;bcb;bct;bcf;bcz};
39+
A = - lapGral3D(k,m,dx,n,dy,o,dz,dc,nc);
40+
b = 2*sin(2*pi*Y).*Z.*(1+2*pi^2*X.*(1-X));
41+
b = reshape(b,[],1);
42+
[A0,b0] = addGralBC3D(A,b,k,m,dx,n,dy,o,dz,dc,nc,v);
43+
ua = A0\b0; % approximate solution
44+
ua = reshape(ua,m+2,n,o+2);
45+
46+
47+
% plot slices as surfaces
48+
figure(bvp)
49+
surf(squeeze(Y((m+3)/2,:,:)),squeeze(Z((m+3)/2,:,:)),squeeze(ua((m+3)/2,:,:)));
50+
title('Approximate Solution: 3D Poisson with Periodic BC (Middle YZ slice)');
51+
xlabel('Y');
52+
ylabel('Z');
53+
shading interp;
54+
figure(bvp+10)
55+
surf(squeeze(Y((m+3)/2,:,:)),squeeze(Z((m+3)/2,:,:)),squeeze(ue((m+3)/2,:,:)));
56+
title('Exact Solution: 3D Poisson with Periodic BC (Middle YZ slice)');
57+
xlabel('Y');
58+
ylabel('Z');
59+
shading interp;
60+
figure(bvp+20)
61+
surf(squeeze(X(:,(n+3)/2,:)),squeeze(Z(:,(n+3)/2,:)),squeeze(ua(:,(n+3)/2,:)));
62+
title('Approximate Solution: 3D Poisson with Periodic BC (Middle XZ slice)');
63+
xlabel('X');
64+
ylabel('Z');
65+
shading interp;
66+
figure(bvp+30)
67+
surf(squeeze(X(:,(n+3)/2,:)),squeeze(Z(:,(n+3)/2,:)),squeeze(ue(:,(n+3)/2,:)));
68+
title('Exact Solution: 3D Poisson with Periodic BC (Middle XZ slice)');
69+
xlabel('X');
70+
ylabel('Z');
71+
shading interp;
72+
figure(bvp+40)
73+
surf(squeeze(X(:,:,(o+3)/2)),squeeze(Y(:,:,(o+3)/2)),squeeze(ua(:,:,(o+3)/2)));
74+
title('Approximate Solution: 3D Poisson with Periodic BC (Middle XY slice)');
75+
xlabel('X');
76+
ylabel('Y');
77+
shading interp;
78+
figure(bvp+50)
79+
surf(squeeze(X(:,:,(o+3)/2)),squeeze(Y(:,:,(o+3)/2)),squeeze(ue(:,:,(o+3)/2)));
80+
title('Exact Solution: 3D Poisson with Periodic BC (Middle XY slice)');
81+
xlabel('X');
82+
ylabel('Y');
83+
shading interp;
84+
85+
fprintf('Maximum error: %.4f\n', max(max(max(abs(ue-ua)))))
86+
fprintf('Relative error: %.4f%%\n', 100*max(max(max(abs(ue-ua))))/(max(max(max(ue))) - min(min(min(ue)))))

‎examples/matlab/elliptic3DXDirichletYPerZPer.m‎

Lines changed: 85 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,85 @@
1+
% ====================== Test 2 =====================
2+
% 3D Poisson BVP: Dirichlet (on left and right faces), Periodic, Periodic BC
3+
% -(u_xx + u_yy + u_zz) = - f(x,y,z), -1 < x < 1, 0 < y,z < 2 pi
4+
% - f(x,y,z) = 16(1-x^2)cos(4x) - 16x sin(4x) + 2(cos(4x)+sin(2y)+sin(4z)) + (1-x^2)(4 sin(2y) + 16 sin(4z))
5+
% BC: u(-1,y,z) = 0, u(1,x,y) = 0, u(x,0,z) = u(x,2 pi,z), u_y(x,0,z) = u_y(x,2 pi,z), u(x,y,0) = u(x,y,2 pi), u_z(x,y,0) = u_z(x,y,2 pi)
6+
% exact solution: (cos(4x)+sin(2y)+sin(4z))*(1-x^2)
7+
% ===================================================
8+
% example that uses addGral3D
9+
%
10+
close all; clc;
11+
12+
addpath('../../src/matlab');
13+
14+
k = 2;
15+
bvp = 2;
16+
m = 29; % should be odd to be able to plot the middle slice
17+
n = m+2; % should be odd to be able to plot the middle slice
18+
o = m+4; % should be odd to be able to plot the middle slice
19+
dx = 2/m;
20+
dy = 2*pi/n;
21+
dz = 2*pi/o;
22+
% centers and vertices
23+
xc = [-1 -1+dx/2:dx:1-dx/2 1]';
24+
yc = (dy/2:dy:2*pi-dy/2)';
25+
zc = (dz/2:dz:2*pi-dz/2)';
26+
[Y,X,Z] = meshgrid(yc,xc,zc);
27+
t = 'u_xx + u_yy + u_zz = f(x,y,z), -1 < x < 1, 0 < y,z < 2 pi, u(-1,y,z) = 0, u(1,y,z) = 0, periodic BC on y,z, with exact solution u(x,y,z) = (cos(4x)+sin(2y)+sin(4z))*(1-x^2)';
28+
ue = (cos(4*X)+sin(2*Y)+sin(4*Z)).*(1-X.^2);
29+
dc = [1;1;0;0;0;0];
30+
nc = [0;0;0;0;0;0];
31+
bcl = zeros(n*o,1);
32+
bcr = zeros(n*o,1);
33+
bcb = 0; % zeros((m+2)*o,1);
34+
bct = 0; % zeros((m+2)*o,1);
35+
bcf = 0; % zeros((n+2)*(m+2),1);
36+
bcz = 0; % zeros((n+2)*(m+2),1);
37+
v = {bcl;bcr;bcb;bct;bcf;bcz};
38+
A = - lapGral3D(k,m,dx,n,dy,o,dz,dc,nc);
39+
b = 16*(1-X.^2).*cos(4*X) - 16*X.*sin(4*X) + 2*(cos(4*X)+sin(2*Y)+sin(4*Z)) + (1-X.^2).*(4*sin(2*Y) + 16*sin(4*Z));
40+
b = reshape(b,[],1);
41+
[A0,b0] = addGralBC3D(A,b,k,m,dx,n,dy,o,dz,dc,nc,v);
42+
ua = A0\b0; % approximate solution
43+
ua = reshape(ua,m+2,n,o);
44+
45+
46+
% plot slices as surfaces
47+
figure(bvp)
48+
surf(squeeze(Y((m+1)/2,:,:)),squeeze(Z((m+1)/2,:,:)),squeeze(ua((m+1)/2,:,:)));
49+
title('Approximate Solution: 3D Poisson with Periodic BC (Middle YZ slice)');
50+
xlabel('Y');
51+
ylabel('Z');
52+
shading interp;
53+
figure(bvp+10)
54+
surf(squeeze(Y((m+1)/2,:,:)),squeeze(Z((m+1)/2,:,:)),squeeze(ue((m+1)/2,:,:)));
55+
title('Exact Solution: 3D Poisson with Periodic BC (Middle YZ slice)');
56+
xlabel('Y');
57+
ylabel('Z');
58+
shading interp;
59+
figure(bvp+20)
60+
surf(squeeze(X(:,(n+1)/2,:)),squeeze(Z(:,(n+1)/2,:)),squeeze(ua(:,(n+1)/2,:)));
61+
title('Approximate Solution: 3D Poisson with Periodic BC (Middle XZ slice)');
62+
xlabel('X');
63+
ylabel('Z');
64+
shading interp;
65+
figure(bvp+30)
66+
surf(squeeze(X(:,(n+1)/2,:)),squeeze(Z(:,(n+1)/2,:)),squeeze(ue(:,(n+1)/2,:)));
67+
title('Exact Solution: 3D Poisson with Periodic BC (Middle XZ slice)');
68+
xlabel('X');
69+
ylabel('Z');
70+
shading interp;
71+
figure(bvp+40)
72+
surf(squeeze(X(:,:,(o+1)/2)),squeeze(Y(:,:,(o+1)/2)),squeeze(ua(:,:,(o+1)/2)));
73+
title('Approximate Solution: 3D Poisson with Periodic BC (Middle XY slice)');
74+
xlabel('X');
75+
ylabel('Y');
76+
shading interp;
77+
figure(bvp+50)
78+
surf(squeeze(X(:,:,(o+1)/2)),squeeze(Y(:,:,(o+1)/2)),squeeze(ue(:,:,(o+1)/2)));
79+
title('Exact Solution: 3D Poisson with Periodic BC (Middle XY slice)');
80+
xlabel('X');
81+
ylabel('Y');
82+
shading interp;
83+
84+
fprintf('Maximum error: %.4f\n', max(max(max(abs(ue-ua)))))
85+
fprintf('Relative error: %.4f%%\n', 100*max(max(max(abs(ue-ua))))/(max(max(max(ue))) - min(min(min(ue)))))

0 commit comments

Comments
 (0)