Skip to content

Commit 0cac0fc

Browse files
mdumettmdumett
committed
Merge branch 'addBCs3' of https://github.com/csrc-sdsu/mole into addBCs3
2 parents d5a3eed + 61dc2d3 commit 0cac0fc

32 files changed

Lines changed: 1591 additions & 0 deletions

‎examples/matlab/elliptic1DHomogeneousDirichlet.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 addBC1D
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 = - lap(k,m,dx);
24+
b = ones(size(A,2),1);
25+
[A0,b0] = addBC1D(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/elliptic1DLeftDirichletRightNeumann.m‎

Lines changed: 34 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,34 @@
1+
% ====================== Test 3 =====================
2+
% 1D Poisson BVP: Dirichlet, Neumann Homogeneous BC
3+
% - u'' = 1, 0 < x < 1, u'(0) = 0, u(1) = 0
4+
% exact solution: u(x) = (1 - x^2)/2
5+
% ===================================================
6+
% example that uses addBC1D
7+
%
8+
close all; clc;
9+
10+
addpath('../../src/matlab');
11+
12+
k = 2;
13+
bvp = 3;
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) = (1-x^2)/2';
19+
ue = 0.5*(1-xc.*xc); % exact solution
20+
dc = [0;1];
21+
nc = [1;0];
22+
v = [0;0];
23+
A = - lap(k,m,dx);
24+
b = ones(size(A,2),1);
25+
[A0,b0] = addBC1D(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/elliptic1DLeftDirichletRightRobin.m‎

Lines changed: 40 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,40 @@
1+
% ====================== Test 7 =====================
2+
% 1D Poisson BVP: Dirichlet, Robin BC
3+
% - u'' = pi^2 sin(pi x), 0 < x < 1, u(0) = c, b u(1)+ u'(1) = d
4+
% exact solution: u(x) = sin(pi x) + Ex + F
5+
% E = (d - bc + pi)/(b+1), F = c
6+
% Taken from
7+
% https://www.scirp.org/journal/paperinformation?paperid=50586
8+
%
9+
% b = 400, c = 10, d = 15
10+
% So, E = (pi - 3985)/401, F = 10
11+
% ===================================================
12+
% example that uses addBC1D
13+
%
14+
close all; clc;
15+
16+
addpath('../../src/matlab');
17+
18+
k = 2;
19+
bvp = 7;
20+
m = 2*k+1;
21+
dx = 1/m;
22+
% centers and vertices
23+
xc = [0 dx/2:dx:1-dx/2 1]';
24+
t = '- u" = sin(pi x), 0 < x < 1, u(0) = 10, 400 u(1) + u''(1) = 15, with exact solution u(x) = sin(pi x) + (pi - 3985)x/401 + 10';
25+
ue = sin(pi*xc) + (pi - 3985)*xc/401 + 10; % exact solution
26+
dc = [1;400];
27+
nc = [0;1];
28+
v = [10;15];
29+
A = - lap(k,m,dx);
30+
b = pi^2 * sin(pi*xc);
31+
[A0,b0] = addBC1D(A,b,k,m,dx,dc,nc,v);
32+
ua = A0\b0; % approximate solution
33+
34+
% plot
35+
figure(bvp)
36+
plot(xc,ue,'b*',xc,ua,'ro');
37+
title(t); %,'interpreter','latex');
38+
xlabel('x');
39+
ylabel('u');
40+
legend({'exact','approx'});

‎examples/matlab/elliptic1DLeftNeumannRightNeumann.m‎

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

‎examples/matlab/elliptic1DLeftNeumannRightRobin.m‎

Lines changed: 40 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,40 @@
1+
% ====================== Test 8 =====================
2+
% 1D Poisson BVP: Neumann, Robin BC
3+
% - u'' = pi^2 sin(pi x), 0 < x < 1, u'(0) = c, b u(1)+ u'(1) = d
4+
% exact solution: u(x) = sin(pi x) + Ex + F
5+
% E = c - pi, F = (d + pi - (b + 1)(c - pi))/b
6+
% Taken from
7+
% https://www.scirp.org/journal/paperinformation?paperid=50586
8+
%
9+
% b = 400, c = 10, d = 15
10+
% So, E = - (10 + pi), F = (402 pi + 4025)/400
11+
% ===================================================
12+
% example that uses addBC1D
13+
%
14+
close all; clc;
15+
16+
addpath('../../src/matlab');
17+
18+
k = 2;
19+
bvp = 8;
20+
m = 2*k+1;
21+
dx = 1/m;
22+
% centers and vertices
23+
xc = [0 dx/2:dx:1-dx/2 1]';
24+
t = '- u^" = sin(pi x), 0 < x < 1, u''(0) = 10, 400 u(1) + u''(1) = 15, with exact solution u(x) = sin(pi x) - (10 + pi)x + (402 pi + 4025)/400';
25+
ue = sin(pi*xc) - (10 + pi)*xc + (402*pi+4025)/400; % exact solution
26+
dc = [0;400];
27+
nc = [1;1];
28+
v = [10;15];
29+
A = - lap(k,m,dx);
30+
b = pi^2 * sin(pi*xc);
31+
[A0,b0] = addBC1D(A,b,k,m,dx,dc,nc,v);
32+
ua = A0\b0; % approximate solution
33+
34+
% plot
35+
figure(bvp)
36+
plot(xc,ue,'b*',xc,ua,'ro');
37+
title(t); %,'interpreter','latex');
38+
xlabel('x');
39+
ylabel('u');
40+
legend({'exact','approx'});

‎examples/matlab/elliptic1DLeftRobinRightRobin.m‎

Lines changed: 42 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,42 @@
1+
% ====================== Test 6 =====================
2+
% 1D Poisson BVP: Robin, Robin BC
3+
% - u'' = pi^2 sin(pi x), 0 < x < 1,
4+
% a u(0) + u'(0) = c, b u(1)+ u'(1) = d
5+
% exact solution: u(x) = sin(pi x) + Ex + F
6+
% E = (bc - ad - (a+b)pi)/(b-a(b+1))
7+
% F = (d - (b+1)c + (b+2)pi)/(b-a(b+1))
8+
% Taken from
9+
% https://www.scirp.org/journal/paperinformation?paperid=50586
10+
%
11+
% a = -200, b = 400, c = 10, d = 15
12+
% So, E = (35 - pi)/403, F = (402 pi - 3995)/80600
13+
% ===================================================
14+
% example that uses addBC1D
15+
%
16+
close all; clc;
17+
18+
addpath('../../src/matlab');
19+
20+
k = 2;
21+
bvp = 6;
22+
m = 20;
23+
dx = 1/m;
24+
% centers and vertices
25+
xc = [0 dx/2:dx:1-dx/2 1]';
26+
t = '- u" = sin(pi x), 0 < x < 1, -200 u(0) + u''(0) = 10, 400 u(1) + u''(1) = 15, with exact solution u(x) = sin(pi x) + (35-pi)x/403 + (402 pi - 3995)/80600';
27+
ue = sin(pi*xc) + (35 - pi)*xc/403 + (402*pi - 3995)/80600; % exact solution
28+
dc = [-200;400];
29+
nc = [1;1];
30+
v = [10;15];
31+
A = - lap(k,m,dx);
32+
b = pi^2 * sin(pi*xc);
33+
[A0,b0] = addBC1D(A,b,k,m,dx,dc,nc,v);
34+
ua = A0\b0; % approximate solution
35+
36+
% plot
37+
figure(bvp)
38+
plot(xc,ue,'b*',xc,ua,'ro');
39+
title(t); %,'interpreter','latex');
40+
xlabel('x');
41+
ylabel('u');
42+
legend({'exact','approx'});

‎examples/matlab/elliptic1DNonHomogeneousDirichlet.m‎

Lines changed: 34 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,34 @@
1+
% ====================== Test 2 =====================
2+
% 1D Poisson BVP: Dirichlet, Dirichlet Non-Homogeneous BC
3+
% - u'' = 1, 0 < x < 1, u(0) = 1/2, u(1) = 1/2
4+
% exact solution: u(x) = (-x^2 + x + 1)/2
5+
% ===================================================
6+
% example that uses addBC1D
7+
%
8+
close all; clc;
9+
10+
addpath('../../src/matlab');
11+
12+
k = 2;
13+
bvp = 2;
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) = 1/2, u(1) = 1/2, with exact solution u(x) = (-x^2 + x +1)/2';
19+
ue = 0.5*(-xc.*xc + xc + 1); % exact solution
20+
dc = [1;1];
21+
nc = [0;0];
22+
v = [1/2;1/2];
23+
A = - lap(k,m,dx);
24+
b = ones(size(A,2),1);
25+
[A0,b0] = addBC1D(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/elliptic1DPeriodicBC.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 addBC1D
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 = [0 dx/2:dx:1-dx/2 1]';
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 = - lap(k,m,dx);
25+
b = 4*pi^2 * sin(2*pi*xc);
26+
[A0,b0] = addBC1D(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/elliptic1DaddBC.m‎

Lines changed: 33 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,33 @@
1+
% Solves the 1D Poisson's equation with Robin boundary conditions
2+
3+
clc
4+
close all
5+
6+
addpath('../../src/matlab')
7+
8+
west = 0; % Domain's limits
9+
east = 1;
10+
11+
k = 6; % Operator's order of accuracy
12+
m = 2*k+1; % Minimum number of cells to attain the desired accuracy
13+
dx = (east-west)/m; % Step length
14+
grid = [west west+dx/2 : dx : east-dx/2 east]; % grid
15+
16+
% Robin boundary conditions
17+
dc = [1;1];
18+
nc = [1;1];
19+
v = [0;2*exp(1)];
20+
21+
L = lap(k, m, dx); % 1D Mimetic laplacian operator
22+
U = exp(grid)'; % RHS
23+
[L0,U0] = addBC1D(L,U,k,m,dx,dc,nc,v); % add BC to system
24+
U0 = L0\U0; % Solve a linear system of equations
25+
26+
% Plot result
27+
plot(grid, U0, 'o')
28+
hold on
29+
plot(grid, exp(grid))
30+
legend('Approximated', 'Analytical', 'Location', 'NorthWest')
31+
title('Poisson''s equation with Robin BC')
32+
xlabel('x')
33+
ylabel('u(x)')

‎examples/matlab/elliptic2DXDirichletYDirichlet.m‎

Lines changed: 51 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,51 @@
1+
% ====================== Test 1 =====================
2+
% 2D Laplace BVP: Dirichlet, Dirichlet
3+
% u_xx + u_yy = 0, 0 < x,y < pi,
4+
% BC: u(x,0) = e^x, u(x,pi) = - e^x, u(0,y) = cos(y), u(pi,y) = e^pi cos(y)
5+
% exact solution: u(x,y) = e^x cos(y)
6+
% ===================================================
7+
% example that uses addBC2D
8+
%
9+
close all; clc;
10+
11+
addpath('../../src/matlab');
12+
13+
k = 2;
14+
bvp = 1;
15+
m = 99; % it should be odd
16+
n = m+2; % it should be odd
17+
dx = pi/m;
18+
dy = pi/n;
19+
% centers and vertices
20+
xc = [0 dx/2:dx:pi-dx/2 pi]';
21+
yc = [0 dy/2:dy:pi-dy/2 pi]';
22+
[Y,X] = meshgrid(yc,xc);
23+
% t = 'u_xx + u_yy = 0, (x,y) in [0,pi]x[0,pi], u(x,0) = e^x, u(x,pi) = - e^x, u(0,y) = cos(y), u(pi,y) = e^pi cos(y), with exact solution u(x,y) = e^x cos(y)';
24+
ue = exp(X).*cos(Y); % exact solution
25+
dc = [1;1;1;1];
26+
nc = [0;0;0;0];
27+
bcl = squeeze(ue(1,:))'; % left bc (y increases)
28+
bcr = squeeze(ue(end,:))'; % right bc (y increases)
29+
bcb = squeeze(ue(:,1)); % bottom bc (x increases)
30+
bct = squeeze(ue(:,end)); % top bc (x increases)
31+
bcl = bcl(2:end-1,1);
32+
bcr = bcr(2:end-1,1);
33+
v = {bcl;bcr;bcb;bct};
34+
A = - lap2D(k,m,dx,n,dy);
35+
b = zeros(m+2,n+2);
36+
b = reshape(b,[],1);
37+
[A0,b0] = addBC2D(A,b,k,m,dx,n,dy,dc,nc,v);
38+
ua = A0\b0; % approximate solution
39+
ua = reshape(ua,m+2,n+2);
40+
41+
figure(bvp)
42+
surf(X,Y,ua);
43+
title('Approximate Solution: 2D Poisson with Periodic BC');
44+
shading interp;
45+
figure(bvp+10)
46+
surf(X,Y,ue);
47+
title('Exact Solution: 2D Poisson with Periodic BC');
48+
shading interp;
49+
50+
fprintf('Maximum error: %.4f\n', max(max(abs(ue-ua))))
51+
fprintf('Relative error: %.4f%%\n', 100*max(max(abs(ue-ua)))/(max(max(ue)) - min(min(ue))))

0 commit comments

Comments
 (0)