Skip to content

Commit a74cb9f

Browse files
mdumettmdumett
committed
updates boundary info into discrete analog system
addBC functionality adds boundary information into the discrete analog of the PDE, both operator and rhs
1 parent f67492e commit a74cb9f

28 files changed

Lines changed: 1351 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/elliptic2DXDirichletYDirichlet.m‎

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

0 commit comments

Comments
 (0)