new examples and name typos

mdumett · mdumett · commit 61dc2d394fe2 · 2025-02-19T18:07:50.000-08:00
Corrected function names and added four examples: two are old ones written with addBC instead of robinBC, one 1D wave equation using addBC and one 3D with Neumann and periodic BC
diff --git a/examples/matlab/elliptic1DaddBC.m b/examples/matlab/elliptic1DaddBC.m
@@ -0,0 +1,33 @@
+% Solves the 1D Poisson's equation with Robin boundary conditions
+
+clc
+close all
+
+addpath('../../src/matlab')
+
+west = 0;  % Domain's limits
+east = 1;
+
+k = 6;  % Operator's order of accuracy
+m = 2*k+1;  % Minimum number of cells to attain the desired accuracy
+dx = (east-west)/m;  % Step length
+grid = [west west+dx/2 : dx : east-dx/2 east]; % grid
+
+% Robin boundary conditions
+dc = [1;1];
+nc = [1;1];
+v = [0;2*exp(1)];
+
+L = lap(k, m, dx);  % 1D Mimetic laplacian operator
+U = exp(grid)'; % RHS
+[L0,U0] = addBC1D(L,U,k,m,dx,dc,nc,v); % add BC to system
+U0 = L0\U0;  % Solve a linear system of equations
+
+% Plot result
+plot(grid, U0, 'o')
+hold on
+plot(grid, exp(grid))
+legend('Approximated', 'Analytical', 'Location', 'NorthWest')
+title('Poisson''s equation with Robin BC')
+xlabel('x')
+ylabel('u(x)')
diff --git a/examples/matlab/elliptic2DXDirichletYDirichlet.m b/examples/matlab/elliptic2DXDirichletYDirichlet.m
@@ -11,7 +11,6 @@
 addpath('../../src/matlab');
 
 k = 2;
-
 bvp = 1;
 m = 99; % it should be odd
 n = m+2; % it should be odd
diff --git a/examples/matlab/elliptic3DXDirichletYDirichletZDirichlet.m b/examples/matlab/elliptic3DXDirichletYDirichletZDirichlet.m
@@ -11,7 +11,6 @@
 addpath('../../src/matlab');
 
 k = 2;
-
 bvp = 1;
 m = 29; % should be odd to be able to plot the middle slice
 n = m+2; % should be odd to be able to plot the middle slice
@@ -50,6 +49,7 @@
 [A0,b0] = addBC3D(A,b,k,m,dx,n,dy,o,dz,dc,nc,v);
 ua = A0\b0; % approximate solution
 ua = reshape(ua,m+2,n+2,o+2);
+ua = ua - ua((m+1)/2,(n+3)/2,(o+3)/2);
 
 % plot slices as surfaces
 figure(bvp)
diff --git a/examples/matlab/elliptic3DXDirichletYPeriodicZPeriodic.m b/examples/matlab/elliptic3DXDirichletYPeriodicZPeriodic.m
@@ -13,7 +13,6 @@
 addpath('../../src/matlab');
 
 k = 2;
-
 bvp = 2;
 m = 49; % should be odd to be able to plot the middle slice
 n = m+2; % should be odd to be able to plot the middle slice
@@ -43,7 +42,6 @@
 [A0,b0] = addBC3D(A,b,k,m,dx,n,dy,o,dz,dc,nc,v);
 ua = A0\b0; % approximate solution
 ua = reshape(ua,m+2,n+2,o+2);
-% addBC3DSurfacePlots(X,Y,Z,m,n,o,ua,ue,bvp);
 
 
 % plot slices as surfaces
diff --git a/examples/matlab/elliptic3DXNeumannYPeriodicZPeriodic.m b/examples/matlab/elliptic3DXNeumannYPeriodicZPeriodic.m
@@ -0,0 +1,92 @@
+% ====================== Test 2 =====================
+% 3D Poisson BVP: Robin (on left and right faces), Periodic, Periodic BC
+% -(u_xx + u_yy + u_zz) = - f(x,y,z), -1 < x < 1, 0 < y,z < 2 pi
+% - 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))
+% % f(x,y,z) = 16(1-x^2)cos(4x) - 16x sin(4x) + 2 cos(4x) + (1-x^2)(4 sin(2y) + 16 sin(4z))
+% 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)
+% exact solution: (cos(4x)+sin(2y)+sin(4z))*(1-x^2)
+% ===================================================
+% example that uses addBC3D
+%
+close all; clc;
+
+addpath('../../src/matlab');
+
+k = 2;
+bvp = 2;
+m = 39; % should be odd to be able to plot the middle slice
+n = m+2; % should be odd to be able to plot the middle slice
+o = m+4; % should be odd to be able to plot the middle slice
+dx = 2/m;
+dy = 2*pi/n;
+dz = 2*pi/o;
+% centers and vertices
+xc = [-1 -1+dx/2:dx:1-dx/2 1]';
+yc = [0 dy/2:dy:2*pi-dy/2 2*pi]';
+zc = [0 dz/2:dz:2*pi-dz/2 2*pi]';
+[Y,X,Z] = meshgrid(yc,xc,zc);
+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)';
+ue = (cos(4*X)+sin(2*Y)+sin(4*Z)).*(1-X.^2);
+dc = [0;0;0;0;0;0];
+nc = [1;1;0;0;0;0];
+bcl = -2*(cos(4)+sin(2*Y)+sin(4*Z));
+bcr = -2*(cos(4)+sin(2*Y)+sin(4*Z));
+bcl = squeeze(bcl(1,:,:)); % left bc (y increase along rows, z increase along cols)
+bcr = squeeze(bcr(end,:,:)); % right bc (y increase along rows, z increase along cols)
+bcl = reshape(bcl(2:end-1,2:end-1),[],1);
+bcr = reshape(bcr(2:end-1,2:end-1),[],1);
+bcb = zeros((m+2)*o,1);
+bct = zeros((m+2)*o,1);
+bcf = zeros((n+2)*(m+2),1);
+bcz = zeros((n+2)*(m+2),1);
+v = {bcl;bcr;bcb;bct;bcf;bcz};
+A = - lap3D(k,m,dx,n,dy,o,dz);
+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));
+b = reshape(b,[],1);
+[A0,b0] = addBC3D(A,b,k,m,dx,n,dy,o,dz,dc,nc,v);
+ua = A0\b0; % approximate solution
+ua = reshape(ua,m+2,n+2,o+2);
+ua = ua - ua(1) + ue(1);
+% ua = ua - ua((m+1)/2,(n+3)/2,(o+3)/2);
+
+
+% plot slices as surfaces
+figure(bvp)
+surf(squeeze(Y((m+3)/2,:,:)),squeeze(Z((m+3)/2,:,:)),squeeze(ua((m+3)/2,:,:)));
+title('Approximate Solution: 3D Poisson with Periodic BC (Middle YZ slice)');
+xlabel('Y');
+ylabel('Z');
+shading interp;
+figure(bvp+10)
+surf(squeeze(Y((m+3)/2,:,:)),squeeze(Z((m+3)/2,:,:)),squeeze(ue((m+3)/2,:,:)));
+title('Exact Solution: 3D Poisson with Periodic BC (Middle YZ slice)');
+xlabel('Y');
+ylabel('Z');
+shading interp;
+figure(bvp+20)
+surf(squeeze(X(:,(n+3)/2,:)),squeeze(Z(:,(n+3)/2,:)),squeeze(ua(:,(n+3)/2,:)));
+title('Approximate Solution: 3D Poisson with Periodic BC (Middle XZ slice)');
+xlabel('X');
+ylabel('Z');
+shading interp;
+figure(bvp+30)
+surf(squeeze(X(:,(n+3)/2,:)),squeeze(Z(:,(n+3)/2,:)),squeeze(ue(:,(n+3)/2,:)));
+title('Exact Solution: 3D Poisson with Periodic BC (Middle XZ slice)');
+xlabel('X');
+ylabel('Z');
+shading interp;
+figure(bvp+40)
+surf(squeeze(X(:,:,(o+3)/2)),squeeze(Y(:,:,(o+3)/2)),squeeze(ua(:,:,(o+3)/2)));
+title('Approximate Solution: 3D Poisson with Periodic BC (Middle XY slice)');
+xlabel('X');
+ylabel('Y');
+shading interp;
+figure(bvp+50)
+surf(squeeze(X(:,:,(o+3)/2)),squeeze(Y(:,:,(o+3)/2)),squeeze(ue(:,:,(o+3)/2)));
+title('Exact Solution: 3D Poisson with Periodic BC (Middle XY slice)');
+xlabel('X');
+ylabel('Y');
+shading interp;
+
+fprintf('Maximum error: %.4f\n', max(max(max(abs(ue-ua)))))
+fprintf('Relative error: %.4f%%\n', 100*max(max(max(abs(ue-ua))))/(max(max(max(ue))) - min(min(min(ue)))))
diff --git a/examples/matlab/elliptic3DaddBC.m b/examples/matlab/elliptic3DaddBC.m
@@ -0,0 +1,41 @@
+% 3D Staggering example using a 3D Mimetic laplacian
+
+clc
+close all
+
+addpath('../../src/matlab')
+
+k = 2; % Order of accuracy
+m = 5; % -> 7
+n = 6; % -> 8
+o = 7; % -> 9
+
+% boundary type
+dc = [1;1;1;1;1;1];
+nc = [0;0;0;0;0;0];
+bcl = zeros(n*o,1);
+bcr = zeros(n*o,1);
+bcb = zeros((m+2)*o,1);
+bct = zeros((m+2)*o,1);
+bcf = 100*ones((n+2)*(m+2),1); % Dirichlet boundary condition
+bcz = zeros((n+2)*(m+2),1);
+v = {bcl;bcr;bcb;bct;bcf;bcz};
+
+L = lap3D(k, m, 1, n, 1, o, 1); % 3D Mimetic laplacian operator
+RHS = zeros(m+2, n+2, o+2);
+RHS = reshape(RHS, [], 1);
+
+[L0, RHS0] = addBC3D(L,RHS,k,m,1,n,1,o,1,dc,nc,v); % add BC to linear system
+SOL = L0\RHS0;
+SOL = reshape(SOL, m+2, n+2, o+2);
+
+p = 2; % Page to be displayed
+
+page = SOL(:, :, p);
+
+imagesc(page)
+title([num2str(p) ' page'])
+xlabel('m')
+ylabel('n')
+set(gca, 'YDir', 'Normal')
+colorbar
diff --git a/examples/matlab/sturmLiouvilleBessel.m b/examples/matlab/sturmLiouvilleBessel.m
@@ -22,7 +22,7 @@
 nc = [0;0];
 v = [0;besselj(3,1)];
 G = grad(k,m,dx);
-I = interpolFacesToStaggered1D(k,m);
+I = interpolFacesToCentersG1D(k,m);
 A = sparse(diag(xc.^2)*lap(k,m,dx) + diag(xc)*I*G + diag(xc.^2 - 9)*speye(m+2,m+2)); % nu = 3
 b = zeros(size(A,2),1);
 [A0,b0] = addBC1D(A,b,k,m,dx,dc,nc,v);
diff --git a/examples/matlab/sturmLiouvilleChebyshev.m b/examples/matlab/sturmLiouvilleChebyshev.m
@@ -21,7 +21,7 @@
 nc = [0;0];
 v = [1;1];
 G = grad(k,m,dx);
-I = interpolFacesToStaggered1D(k,m);
+I = interpolFacesToCentersG1D(k,m);
 A = sparse(diag(1-xc.^2)*lap(k,m,dx) - diag(xc)*I*G) + 4*speye(m+2,m+2); % n = 2
 b = zeros(size(A,2),1);
 [A0,b0] = addBC1D(A,b,k,m,dx,dc,nc,v);
diff --git a/examples/matlab/sturmLiouvilleHermite.m b/examples/matlab/sturmLiouvilleHermite.m
@@ -24,7 +24,7 @@
 % v = [hermiteH(4,-1);hermiteH(4,1)];
 v = [-20; -20];
 G = grad(k,m,dx);
-I = interpolFacesToStaggered1D(k,m);
+I = interpolFacesToCentersG1D(k,m);
 A = lap(k,m,dx) - 2*sparse(diag(xc)*I*G) + 8*speye(m+2,m+2); % m = 4
 b = zeros(size(A,2),1);
 [A0,b0] = addBC1D(A,b,k,m,dx,dc,nc,v);
diff --git a/examples/matlab/sturmLiouvilleLaguerre.m b/examples/matlab/sturmLiouvilleLaguerre.m
@@ -30,7 +30,7 @@
 % v = [laguerreL(4,0);laguerreL(4,2)];
 v = [1; 0.333333333333333];
 G = grad(k,m,dx);
-I = interpolFacesToStaggered1D(k,m);
+I = interpolFacesToCentersG1D(k,m);
 A = sparse(diag(xc)*lap(k,m,dx) + diag(1-xc)*I*G) + 4*speye(m+2,m+2); % n = 4
 b = zeros(size(A,2),1);
 [A0,b0] = addBC1D(A,b,k,m,dx,dc,nc,v);
diff --git a/examples/matlab/sturmLiouvilleLegendre.m b/examples/matlab/sturmLiouvilleLegendre.m
@@ -24,7 +24,7 @@
 nc = [0;0];
 v = [-1;1];
 G = grad(k,m,dx);
-I = interpolFacesToStaggered1D(k,m);
+I = interpolFacesToCentersG1D(k,m);
 A = sparse(diag(1-xc.^2)*lap(k,m,dx) - 2*diag(xc)*I*G) + 12*speye(m+2,m+2); % n = 3
 b = zeros(size(A,2),1);
 [A0,b0] = addBC1D(A,b,k,m,dx,dc,nc,v);
diff --git a/examples/matlab/wave1DTimeVaryingBC.m b/examples/matlab/wave1DTimeVaryingBC.m
@@ -0,0 +1,70 @@
+% ====================== Test 2 =====================
+% 1D wave equation 
+% using addBC time central difference and implicit scheme: BC Homogeneous Dirichlet
+% u_tt = c^2 u_xx, 0 < x < M, 0 < t < T, 
+% initial condition: u(0,t) = sin(pi x), u_t(0,t) = 0
+% boundary condition: u(0,t) = 0, u(M,t) = 0
+% exact solution: sin(pi x) cos(2 pi t)
+%
+% Taken from example 1 of 
+% https://www.uni-muenster.de/Physik.TP/archive/fileadmin/lehre/NumMethoden/WS1011/script1011Wave.pdf
+% ===================================================
+% example that uses addBC1D
+%
+% scheme:
+% ((U^(n+1) - 2 U^(n) + U^(n-1))/dt = (c^2/2) L (U^(n+1) + U^(n-1))
+%
+clear all; close all; clc
+
+addpath('../../src/matlab')
+
+% some other parameters
+c = 2;
+fx = 10; % final x
+ft = 2; % final time
+
+% Spatial discretization
+k = 2;     % Order of accuracy (spatial) 
+bvp = 9;   % figure index
+m = 300;   % Number of cells
+dx = fx/m;  % Step length
+dt = dx/c; % Time length (c dt/dx =< 1 or dt =< dx/c)
+
+% centers and vertices
+xc = [0 dx/2:dx:fx-dx/2 fx]';
+ts = (0:dt:ft)';
+[T,X] = meshgrid(ts,xc);
+t = 'u_tt = 4 u_xx, 0 < x < M, 0 < t < T, homogeneous Dirichlet BC, u(0,t) = sin(pi x), u_t(0,t) = 0, with exact solution u(x) = sin(pi x) cos(2 pi t)';
+ue = sin(pi*X).*cos(2*pi*T); % exact solution
+% boundary condition does not depend on time, otherwise should go inside time loop
+dc = [1;1];
+nc = [0;0];
+v = [0;0];
+% implicit scheme
+a = (dt*c)^2/2;
+% (I - a L)U^(n+1) = 2 U^n - (I - a L)U^(n-1)
+A = speye(m+2,m+2) - a*lap(k,m,dx);
+% approximate solution storage
+ua = zeros(size(ue));
+ua(:,1) = sin(pi*xc); % initial position
+ua(:,2) = ue(:,2); % from exact solution instead of using initial velocity
+% time loop
+for idx = 3:ft/dt + 1
+    b = 2*ua(:,idx-1) - A*ua(:,idx-2);
+    [A0,b0] = addBC1D(A,b,k,m,dx,dc,nc,v);
+    ua(:,idx) = A0\b0; % approximate solution
+    % ua0 = ua(:,idx)'
+    % ue0 = ue(:,idx)'
+end
+
+figure(bvp)
+surf(X,T,ua);
+title('Approximate Solution: 1D Wave Equation with Homogeneous Dirichlet BC');
+shading interp;
+figure(bvp+10)
+surf(X,T,ue);
+title('Exact Solution: 1D Wave Equation with Homogeneous Dirichlet BC');
+shading interp;
+
+fprintf('Maximum error: %.4f\n', max(max(abs(ue-ua))))
+fprintf('Relative error: %.4f%%\n', 100*max(max(abs(ue-ua)))/(max(max(ue)) - min(min(ue))))
diff --git a/src/matlab/addBC3D.m b/src/matlab/addBC3D.m
@@ -87,16 +87,23 @@
     [Abcl,Abcr,Abcb,Abct,Abcf,Abcz] = addBC3Dlhs(k, m, dx, n, dy, o, dz, dc, nc);
 
     % get rhs entries affected by bcs for left, right, bottom, top, front, back faces, resp.
-    [rl,~,~] = find(Abcl);
-    [rr,~,~] = find(Abcr);
-    [rb,~,~] = find(Abcb);
-    [rt,~,~] = find(Abct);
-    [rf,~,~] = find(Abcf);
-    [rz,~,~] = find(Abcz);
+    [rl,~,~] = find(Abcl); 
+    [rr,~,~] = find(Abcr); 
+    [rb,~,~] = find(Abcb); 
+    [rt,~,~] = find(Abct); 
+    [rf,~,~] = find(Abcf); 
+    [rz,~,~] = find(Abcz); 
+    rl = unique(rl);
+    rr = unique(rr);
+    rb = unique(rb);
+    rt = unique(rt);
+    rf = unique(rf);
+    rz = unique(rz);
 
     % remove rows from matrix A
     Abc = Abcl + Abcr + Abcb + Abct + Abcz + Abcf;
     [rowsbc,~,~] = find(Abc);
+    rowsbc = unique(rowsbc);
     [rows,cols,s] = find(A(rowsbc,:));
     A = A - sparse(rows, cols, s, size(A,1), size(A,2));
     % update matrix A with boundary information
