Mercurial > repos > public > sbplib
changeset 1286:4cb627c7fb90 feature/boundary_optimized_grids
Make D1Nonequidistant use the grid generation functions accurate/minimalBoundaryOptimizedGrid and remove grid generation from +implementations
| author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
|---|---|
| date | Wed, 01 Jul 2020 13:43:32 +0200 |
| parents | 6b68f939d023 |
| children | 38653d26225c |
| files | +sbp/+implementations/d1_noneq_10.m +sbp/+implementations/d1_noneq_12.m +sbp/+implementations/d1_noneq_4.m +sbp/+implementations/d1_noneq_6.m +sbp/+implementations/d1_noneq_8.m +sbp/+implementations/d1_noneq_minimal_10.m +sbp/+implementations/d1_noneq_minimal_12.m +sbp/+implementations/d1_noneq_minimal_4.m +sbp/+implementations/d1_noneq_minimal_6.m +sbp/+implementations/d1_noneq_minimal_8.m +sbp/D1Nonequidistant.m |
| diffstat | 11 files changed, 54 insertions(+), 516 deletions(-) [+] |
line wrap: on
line diff
--- a/+sbp/+implementations/d1_noneq_10.m Wed Jul 01 11:15:57 2020 +0200 +++ b/+sbp/+implementations/d1_noneq_10.m Wed Jul 01 13:43:32 2020 +0200 @@ -1,48 +1,12 @@ -function [D1,H,x,h] = d1_noneq_10(N,L) +function [D1,H] = d1_noneq_10(N,h) -% L: Domain length % N: Number of grid points -if(nargin < 2) - L = 1; -end - if(N<20) error('Operator requires at least 20 grid points'); end % BP: Number of boundary points -% m: Number of nonequidistant spacings -% order: Accuracy of interior stencil BP = 10; -m = 5; -order = 10; - -%%%% Non-equidistant grid points %%%%% -x0 = 0.0000000000000e+00; -x1 = 3.5902433622052e-01; -x2 = 1.1436659188355e+00; -x3 = 2.2144895894456e+00; -x4 = 3.3682742337736e+00; -x5 = 4.4309689056870e+00; -x6 = 5.4309689056870e+00; -x7 = 6.4309689056870e+00; -x8 = 7.4309689056870e+00; -x9 = 8.4309689056870e+00; -x10 = 9.4309689056870e+00; - -xb = sparse(m+1,1); -for i = 0:m - xb(i+1) = eval(['x' num2str(i)]); -end -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -%%%% Compute h %%%%%%%%%% -h = L/(2*xb(end) + N-1-2*m); -%%%%%%%%%%%%%%%%%%%%%%%%% - -%%%% Define grid %%%%%%%% -x = h*[xb; linspace(xb(end)+1,L/h-xb(end)-1,N-2*(m+1))'; L/h-flip(xb) ]; -%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Norm matrix %%%%%%%% P = sparse(BP,1); @@ -69,22 +33,9 @@ %%%%%%%%%%%%%%%%%%%%%%%%% %%%% Q matrix %%%%%%%%%%% - % interior stencil -switch order - case 2 - d = [-1/2,0,1/2]; - case 4 - d = [1/12,-2/3,0,2/3,-1/12]; - case 6 - d = [-1/60,3/20,-3/4,0,3/4,-3/20,1/60]; - case 8 - d = [1/280,-4/105,1/5,-4/5,0,4/5,-1/5,4/105,-1/280]; - case 10 - d = [-1/1260,5/504,-5/84,5/21,-5/6,0,5/6,-5/21,5/84,-5/504,1/1260]; - case 12 - d = [1/5544,-1/385,1/56,-5/63,15/56,-6/7,0,6/7,-15/56,5/63,-1/56,1/385,-1/5544]; -end +order = 10; +d = [-1/1260,5/504,-5/84,5/21,-5/6,0,5/6,-5/21,5/84,-5/504,1/1260]; d = repmat(d,N,1); Q = spdiags(d,-order/2:order/2,N,N);
--- a/+sbp/+implementations/d1_noneq_12.m Wed Jul 01 11:15:57 2020 +0200 +++ b/+sbp/+implementations/d1_noneq_12.m Wed Jul 01 13:43:32 2020 +0200 @@ -1,50 +1,12 @@ -function [D1,H,x,h] = d1_noneq_12(N,L) +function [D1,H] = d1_noneq_12(N,h) -% L: Domain length % N: Number of grid points -if(nargin < 2) - L = 1; -end - if(N<24) error('Operator requires at least 24 grid points'); end % BP: Number of boundary points -% m: Number of nonequidistant spacings -% order: Accuracy of interior stencil BP = 12; -m = 6; -order = 12; - -%%%% Non-equidistant grid points %%%%% -x0 = 0.0000000000000e+00; -x1 = 3.6098032343909e-01; -x2 = 1.1634317168086e+00; -x3 = 2.2975905356987e+00; -x4 = 3.6057529790929e+00; -x5 = 4.8918275675510e+00; -x6 = 6.0000000000000e+00; -x7 = 7.0000000000000e+00; -x8 = 8.0000000000000e+00; -x9 = 9.0000000000000e+00; -x10 = 1.0000000000000e+01; -x11 = 1.1000000000000e+01; -x12 = 1.2000000000000e+01; - -xb = sparse(m+1,1); -for i = 0:m - xb(i+1) = eval(['x' num2str(i)]); -end -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -%%%% Compute h %%%%%%%%%% -h = L/(2*xb(end) + N-1-2*m); -%%%%%%%%%%%%%%%%%%%%%%%%% - -%%%% Define grid %%%%%%%% -x = h*[xb; linspace(xb(end)+1,L/h-xb(end)-1,N-2*(m+1))'; L/h-flip(xb) ]; -%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Norm matrix %%%%%%%% P = sparse(BP,1); @@ -73,22 +35,9 @@ %%%%%%%%%%%%%%%%%%%%%%%%% %%%% Q matrix %%%%%%%%%%% - % interior stencil -switch order - case 2 - d = [-1/2,0,1/2]; - case 4 - d = [1/12,-2/3,0,2/3,-1/12]; - case 6 - d = [-1/60,3/20,-3/4,0,3/4,-3/20,1/60]; - case 8 - d = [1/280,-4/105,1/5,-4/5,0,4/5,-1/5,4/105,-1/280]; - case 10 - d = [-1/1260,5/504,-5/84,5/21,-5/6,0,5/6,-5/21,5/84,-5/504,1/1260]; - case 12 - d = [1/5544,-1/385,1/56,-5/63,15/56,-6/7,0,6/7,-15/56,5/63,-1/56,1/385,-1/5544]; -end +order = 12; +d = [1/5544,-1/385,1/56,-5/63,15/56,-6/7,0,6/7,-15/56,5/63,-1/56,1/385,-1/5544]; d = repmat(d,N,1); Q = spdiags(d,-order/2:order/2,N,N);
--- a/+sbp/+implementations/d1_noneq_4.m Wed Jul 01 11:15:57 2020 +0200 +++ b/+sbp/+implementations/d1_noneq_4.m Wed Jul 01 13:43:32 2020 +0200 @@ -1,42 +1,12 @@ -function [D1,H,x,h] = d1_noneq_4(N,L) +function [D1,H] = d1_noneq_4(N,h) -% L: Domain length % N: Number of grid points -if(nargin < 2) - L = 1; -end - if(N<8) error('Operator requires at least 8 grid points'); end % BP: Number of boundary points -% m: Number of nonequidistant spacings -% order: Accuracy of interior stencil BP = 4; -m = 2; -order = 4; - -%%%% Non-equidistant grid points %%%%% -x0 = 0.0000000000000e+00; -x1 = 6.8764546205559e-01; -x2 = 1.8022115125776e+00; -x3 = 2.8022115125776e+00; -x4 = 3.8022115125776e+00; - -xb = sparse(m+1,1); -for i = 0:m - xb(i+1) = eval(['x' num2str(i)]); -end -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -%%%% Compute h %%%%%%%%%% -h = L/(2*xb(end) + N-1-2*m); -%%%%%%%%%%%%%%%%%%%%%%%%% - -%%%% Define grid %%%%%%%% -x = h*[xb; linspace(xb(end)+1,L/h-xb(end)-1,N-2*(m+1))'; L/h-flip(xb) ]; -%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Norm matrix %%%%%%%% P = sparse(BP,1); @@ -57,22 +27,9 @@ %%%%%%%%%%%%%%%%%%%%%%%%% %%%% Q matrix %%%%%%%%%%% - % interior stencil -switch order - case 2 - d = [-1/2,0,1/2]; - case 4 - d = [1/12,-2/3,0,2/3,-1/12]; - case 6 - d = [-1/60,3/20,-3/4,0,3/4,-3/20,1/60]; - case 8 - d = [1/280,-4/105,1/5,-4/5,0,4/5,-1/5,4/105,-1/280]; - case 10 - d = [-1/1260,5/504,-5/84,5/21,-5/6,0,5/6,-5/21,5/84,-5/504,1/1260]; - case 12 - d = [1/5544,-1/385,1/56,-5/63,15/56,-6/7,0,6/7,-15/56,5/63,-1/56,1/385,-1/5544]; -end +order = 4; +d = [1/12,-2/3,0,2/3,-1/12]; d = repmat(d,N,1); Q = spdiags(d,-order/2:order/2,N,N);
--- a/+sbp/+implementations/d1_noneq_6.m Wed Jul 01 11:15:57 2020 +0200 +++ b/+sbp/+implementations/d1_noneq_6.m Wed Jul 01 13:43:32 2020 +0200 @@ -1,44 +1,12 @@ -function [D1,H,x,h] = d1_noneq_6(N,L) +function [D1,H] = d1_noneq_6(N,h) -% L: Domain length % N: Number of grid points -if(nargin < 2) - L = 1; -end - if(N<12) error('Operator requires at least 12 grid points'); end % BP: Number of boundary points -% m: Number of nonequidistant spacings -% order: Accuracy of interior stencil BP = 6; -m = 3; -order = 6; - -%%%% Non-equidistant grid points %%%%% -x0 = 0.0000000000000e+00; -x1 = 4.4090263368623e-01; -x2 = 1.2855984345073e+00; -x3 = 2.2638953951239e+00; -x4 = 3.2638953951239e+00; -x5 = 4.2638953951239e+00; -x6 = 5.2638953951239e+00; - -xb = sparse(m+1,1); -for i = 0:m - xb(i+1) = eval(['x' num2str(i)]); -end -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -%%%% Compute h %%%%%%%%%% -h = L/(2*xb(end) + N-1-2*m); -%%%%%%%%%%%%%%%%%%%%%%%%% - -%%%% Define grid %%%%%%%% -x = h*[xb; linspace(xb(end)+1,L/h-xb(end)-1,N-2*(m+1))'; L/h-flip(xb) ]; -%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Norm matrix %%%%%%%% P = sparse(BP,1); @@ -61,22 +29,9 @@ %%%%%%%%%%%%%%%%%%%%%%%%% %%%% Q matrix %%%%%%%%%%% - % interior stencil -switch order - case 2 - d = [-1/2,0,1/2]; - case 4 - d = [1/12,-2/3,0,2/3,-1/12]; - case 6 - d = [-1/60,3/20,-3/4,0,3/4,-3/20,1/60]; - case 8 - d = [1/280,-4/105,1/5,-4/5,0,4/5,-1/5,4/105,-1/280]; - case 10 - d = [-1/1260,5/504,-5/84,5/21,-5/6,0,5/6,-5/21,5/84,-5/504,1/1260]; - case 12 - d = [1/5544,-1/385,1/56,-5/63,15/56,-6/7,0,6/7,-15/56,5/63,-1/56,1/385,-1/5544]; -end +order = 6; +d = [-1/60,3/20,-3/4,0,3/4,-3/20,1/60]; d = repmat(d,N,1); Q = spdiags(d,-order/2:order/2,N,N);
--- a/+sbp/+implementations/d1_noneq_8.m Wed Jul 01 11:15:57 2020 +0200 +++ b/+sbp/+implementations/d1_noneq_8.m Wed Jul 01 13:43:32 2020 +0200 @@ -1,46 +1,12 @@ -function [D1,H,x,h] = d1_noneq_8(N,L) +function [D1,H] = d1_noneq_8(N,h) -% L: Domain length % N: Number of grid points -if(nargin < 2) - L = 1; -end - if(N<16) error('Operator requires at least 16 grid points'); end % BP: Number of boundary points -% m: Number of nonequidistant spacings -% order: Accuracy of interior stencil BP = 8; -m = 4; -order = 8; - -%%%% Non-equidistant grid points %%%%% -x0 = 0.0000000000000e+00; -x1 = 3.8118550247622e-01; -x2 = 1.1899550868338e+00; -x3 = 2.2476300175641e+00; -x4 = 3.3192851303204e+00; -x5 = 4.3192851303204e+00; -x6 = 5.3192851303204e+00; -x7 = 6.3192851303204e+00; -x8 = 7.3192851303204e+00; - -xb = sparse(m+1,1); -for i = 0:m - xb(i+1) = eval(['x' num2str(i)]); -end -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -%%%% Compute h %%%%%%%%%% -h = L/(2*xb(end) + N-1-2*m); -%%%%%%%%%%%%%%%%%%%%%%%%% - -%%%% Define grid %%%%%%%% -x = h*[xb; linspace(xb(end)+1,L/h-xb(end)-1,N-2*(m+1))'; L/h-flip(xb) ]; -%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Norm matrix %%%%%%%% P = sparse(BP,1); @@ -65,22 +31,9 @@ %%%%%%%%%%%%%%%%%%%%%%%%% %%%% Q matrix %%%%%%%%%%% - % interior stencil -switch order - case 2 - d = [-1/2,0,1/2]; - case 4 - d = [1/12,-2/3,0,2/3,-1/12]; - case 6 - d = [-1/60,3/20,-3/4,0,3/4,-3/20,1/60]; - case 8 - d = [1/280,-4/105,1/5,-4/5,0,4/5,-1/5,4/105,-1/280]; - case 10 - d = [-1/1260,5/504,-5/84,5/21,-5/6,0,5/6,-5/21,5/84,-5/504,1/1260]; - case 12 - d = [1/5544,-1/385,1/56,-5/63,15/56,-6/7,0,6/7,-15/56,5/63,-1/56,1/385,-1/5544]; -end +order = 8; +d = [1/280,-4/105,1/5,-4/5,0,4/5,-1/5,4/105,-1/280]; d = repmat(d,N,1); Q = spdiags(d,-order/2:order/2,N,N);
--- a/+sbp/+implementations/d1_noneq_minimal_10.m Wed Jul 01 11:15:57 2020 +0200 +++ b/+sbp/+implementations/d1_noneq_minimal_10.m Wed Jul 01 13:43:32 2020 +0200 @@ -1,46 +1,12 @@ -function [D1,H,x,h] = d1_noneq_minimal_10(N,L) +function [D1,H] = d1_noneq_minimal_10(N,h) -% L: Domain length % N: Number of grid points -if(nargin < 2) - L = 1; -end - if(N<16) error('Operator requires at least 16 grid points'); end % BP: Number of boundary points -% m: Number of nonequidistant spacings -% order: Accuracy of interior stencil BP = 8; -m = 3; -order = 10; - -%%%% Non-equidistant grid points %%%%% -x0 = 0.0000000000000e+00; -x1 = 5.8556160757529e-01; -x2 = 1.7473267488572e+00; -x3 = 3.0000000000000e+00; -x4 = 4.0000000000000e+00; -x5 = 5.0000000000000e+00; -x6 = 6.0000000000000e+00; -x7 = 7.0000000000000e+00; -x8 = 8.0000000000000e+00; - -xb = sparse(m+1,1); -for i = 0:m - xb(i+1) = eval(['x' num2str(i)]); -end -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -%%%% Compute h %%%%%%%%%% -h = L/(2*xb(end) + N-1-2*m); -%%%%%%%%%%%%%%%%%%%%%%%%% - -%%%% Define grid %%%%%%%% -x = h*[xb; linspace(xb(end)+1,L/h-xb(end)-1,N-2*(m+1))'; L/h-flip(xb) ]; -%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Norm matrix %%%%%%%% P = sparse(BP,1); @@ -65,22 +31,9 @@ %%%%%%%%%%%%%%%%%%%%%%%%% %%%% Q matrix %%%%%%%%%%% - % interior stencil -switch order - case 2 - d = [-1/2,0,1/2]; - case 4 - d = [1/12,-2/3,0,2/3,-1/12]; - case 6 - d = [-1/60,3/20,-3/4,0,3/4,-3/20,1/60]; - case 8 - d = [1/280,-4/105,1/5,-4/5,0,4/5,-1/5,4/105,-1/280]; - case 10 - d = [-1/1260,5/504,-5/84,5/21,-5/6,0,5/6,-5/21,5/84,-5/504,1/1260]; - case 12 - d = [1/5544,-1/385,1/56,-5/63,15/56,-6/7,0,6/7,-15/56,5/63,-1/56,1/385,-1/5544]; -end +order = 10; +d = [-1/1260,5/504,-5/84,5/21,-5/6,0,5/6,-5/21,5/84,-5/504,1/1260]; d = repmat(d,N,1); Q = spdiags(d,-order/2:order/2,N,N);
--- a/+sbp/+implementations/d1_noneq_minimal_12.m Wed Jul 01 11:15:57 2020 +0200 +++ b/+sbp/+implementations/d1_noneq_minimal_12.m Wed Jul 01 13:43:32 2020 +0200 @@ -1,48 +1,12 @@ -function [D1,H,x,h] = d1_noneq_minimal_12(N,L) +function [D1,H] = d1_noneq_minimal_12(N,h) -% L: Domain length % N: Number of grid points -if(nargin < 2) - L = 1; -end - if(N<20) error('Operator requires at least 20 grid points'); end % BP: Number of boundary points -% m: Number of nonequidistant spacings -% order: Accuracy of interior stencil BP = 10; -m = 4; -order = 12; - -%%%% Non-equidistant grid points %%%%% -x0 = 0.0000000000000e+00; -x1 = 4.6552112904489e-01; -x2 = 1.4647984306493e+00; -x3 = 2.7620429464763e+00; -x4 = 4.0000000000000e+00; -x5 = 5.0000000000000e+00; -x6 = 6.0000000000000e+00; -x7 = 7.0000000000000e+00; -x8 = 8.0000000000000e+00; -x9 = 9.0000000000000e+00; -x10 = 1.0000000000000e+01; - -xb = sparse(m+1,1); -for i = 0:m - xb(i+1) = eval(['x' num2str(i)]); -end -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -%%%% Compute h %%%%%%%%%% -h = L/(2*xb(end) + N-1-2*m); -%%%%%%%%%%%%%%%%%%%%%%%%% - -%%%% Define grid %%%%%%%% -x = h*[xb; linspace(xb(end)+1,L/h-xb(end)-1,N-2*(m+1))'; L/h-flip(xb) ]; -%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Norm matrix %%%%%%%% P = sparse(BP,1); @@ -69,22 +33,9 @@ %%%%%%%%%%%%%%%%%%%%%%%%% %%%% Q matrix %%%%%%%%%%% - % interior stencil -switch order - case 2 - d = [-1/2,0,1/2]; - case 4 - d = [1/12,-2/3,0,2/3,-1/12]; - case 6 - d = [-1/60,3/20,-3/4,0,3/4,-3/20,1/60]; - case 8 - d = [1/280,-4/105,1/5,-4/5,0,4/5,-1/5,4/105,-1/280]; - case 10 - d = [-1/1260,5/504,-5/84,5/21,-5/6,0,5/6,-5/21,5/84,-5/504,1/1260]; - case 12 - d = [1/5544,-1/385,1/56,-5/63,15/56,-6/7,0,6/7,-15/56,5/63,-1/56,1/385,-1/5544]; -end +order = 12; +d = [1/5544,-1/385,1/56,-5/63,15/56,-6/7,0,6/7,-15/56,5/63,-1/56,1/385,-1/5544]; d = repmat(d,N,1); Q = spdiags(d,-order/2:order/2,N,N);
--- a/+sbp/+implementations/d1_noneq_minimal_4.m Wed Jul 01 11:15:57 2020 +0200 +++ b/+sbp/+implementations/d1_noneq_minimal_4.m Wed Jul 01 13:43:32 2020 +0200 @@ -1,41 +1,12 @@ -function [D1,H,x,h] = d1_noneq_minimal_4(N,L) +function [D1,H] = d1_noneq_minimal_4(N,h) -% L: Domain length % N: Number of grid points -if(nargin < 2) - L = 1; -end - if(N<6) error('Operator requires at least 6 grid points'); end % BP: Number of boundary points -% m: Number of nonequidistant spacings -% order: Accuracy of interior stencil BP = 3; -m = 1; -order = 4; - -%%%% Non-equidistant grid points %%%%% -x0 = 0.0000000000000e+00; -x1 = 7.7122987842562e-01; -x2 = 1.7712298784256e+00; -x3 = 2.7712298784256e+00; - -xb = sparse(m+1,1); -for i = 0:m - xb(i+1) = eval(['x' num2str(i)]); -end -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -%%%% Compute h %%%%%%%%%% -h = L/(2*xb(end) + N-1-2*m); -%%%%%%%%%%%%%%%%%%%%%%%%% - -%%%% Define grid %%%%%%%% -x = h*[xb; linspace(xb(end)+1,L/h-xb(end)-1,N-2*(m+1))'; L/h-flip(xb) ]; -%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Norm matrix %%%%%%%% P = sparse(BP,1); @@ -55,22 +26,9 @@ %%%%%%%%%%%%%%%%%%%%%%%%% %%%% Q matrix %%%%%%%%%%% - % interior stencil -switch order - case 2 - d = [-1/2,0,1/2]; - case 4 - d = [1/12,-2/3,0,2/3,-1/12]; - case 6 - d = [-1/60,3/20,-3/4,0,3/4,-3/20,1/60]; - case 8 - d = [1/280,-4/105,1/5,-4/5,0,4/5,-1/5,4/105,-1/280]; - case 10 - d = [-1/1260,5/504,-5/84,5/21,-5/6,0,5/6,-5/21,5/84,-5/504,1/1260]; - case 12 - d = [1/5544,-1/385,1/56,-5/63,15/56,-6/7,0,6/7,-15/56,5/63,-1/56,1/385,-1/5544]; -end +order = 4; +d = [1/12,-2/3,0,2/3,-1/12]; d = repmat(d,N,1); Q = spdiags(d,-order/2:order/2,N,N);
--- a/+sbp/+implementations/d1_noneq_minimal_6.m Wed Jul 01 11:15:57 2020 +0200 +++ b/+sbp/+implementations/d1_noneq_minimal_6.m Wed Jul 01 13:43:32 2020 +0200 @@ -1,43 +1,12 @@ -function [D1,H,x,h] = d1_noneq_minimal_6(N,L) +function [D1,H] = d1_noneq_minimal_6(N,h) -% L: Domain length % N: Number of grid points -if(nargin < 2) - L = 1; -end - if(N<10) error('Operator requires at least 10 grid points'); end % BP: Number of boundary points -% m: Number of nonequidistant spacings -% order: Accuracy of interior stencil BP = 5; -m = 2; -order = 6; - -%%%% Non-equidistant grid points %%%%% -x0 = 0.0000000000000e+00; -x1 = 4.0842950991998e-01; -x2 = 1.1968523189207e+00; -x3 = 2.1968523189207e+00; -x4 = 3.1968523189207e+00; -x5 = 4.1968523189207e+00; - -xb = sparse(m+1,1); -for i = 0:m - xb(i+1) = eval(['x' num2str(i)]); -end -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -%%%% Compute h %%%%%%%%%% -h = L/(2*xb(end) + N-1-2*m); -%%%%%%%%%%%%%%%%%%%%%%%%% - -%%%% Define grid %%%%%%%% -x = h*[xb; linspace(xb(end)+1,L/h-xb(end)-1,N-2*(m+1))'; L/h-flip(xb) ]; -%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Norm matrix %%%%%%%% P = sparse(BP,1); @@ -59,22 +28,9 @@ %%%%%%%%%%%%%%%%%%%%%%%%% %%%% Q matrix %%%%%%%%%%% - % interior stencil -switch order - case 2 - d = [-1/2,0,1/2]; - case 4 - d = [1/12,-2/3,0,2/3,-1/12]; - case 6 - d = [-1/60,3/20,-3/4,0,3/4,-3/20,1/60]; - case 8 - d = [1/280,-4/105,1/5,-4/5,0,4/5,-1/5,4/105,-1/280]; - case 10 - d = [-1/1260,5/504,-5/84,5/21,-5/6,0,5/6,-5/21,5/84,-5/504,1/1260]; - case 12 - d = [1/5544,-1/385,1/56,-5/63,15/56,-6/7,0,6/7,-15/56,5/63,-1/56,1/385,-1/5544]; -end +order = 6; +d = [-1/60,3/20,-3/4,0,3/4,-3/20,1/60]; d = repmat(d,N,1); Q = spdiags(d,-order/2:order/2,N,N);
--- a/+sbp/+implementations/d1_noneq_minimal_8.m Wed Jul 01 11:15:57 2020 +0200 +++ b/+sbp/+implementations/d1_noneq_minimal_8.m Wed Jul 01 13:43:32 2020 +0200 @@ -1,44 +1,12 @@ -function [D1,H,x,h] = d1_noneq_minimal_8(N,L) +function [D1,H] = d1_noneq_minimal_8(N,h) -% L: Domain length % N: Number of grid points -if(nargin < 2) - L = 1; -end - if(N<12) error('Operator requires at least 12 grid points'); end % BP: Number of boundary points -% m: Number of nonequidistant spacings -% order: Accuracy of interior stencil BP = 6; -m = 2; -order = 8; - -%%%% Non-equidistant grid points %%%%% -x0 = 0.0000000000000e+00; -x1 = 4.9439570885261e-01; -x2 = 1.4051531374839e+00; -x3 = 2.4051531374839e+00; -x4 = 3.4051531374839e+00; -x5 = 4.4051531374839e+00; -x6 = 5.4051531374839e+00; - -xb = sparse(m+1,1); -for i = 0:m - xb(i+1) = eval(['x' num2str(i)]); -end -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -%%%% Compute h %%%%%%%%%% -h = L/(2*xb(end) + N-1-2*m); -%%%%%%%%%%%%%%%%%%%%%%%%% - -%%%% Define grid %%%%%%%% -x = h*[xb; linspace(xb(end)+1,L/h-xb(end)-1,N-2*(m+1))'; L/h-flip(xb) ]; -%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Norm matrix %%%%%%%% P = sparse(BP,1); @@ -61,22 +29,9 @@ %%%%%%%%%%%%%%%%%%%%%%%%% %%%% Q matrix %%%%%%%%%%% - % interior stencil -switch order - case 2 - d = [-1/2,0,1/2]; - case 4 - d = [1/12,-2/3,0,2/3,-1/12]; - case 6 - d = [-1/60,3/20,-3/4,0,3/4,-3/20,1/60]; - case 8 - d = [1/280,-4/105,1/5,-4/5,0,4/5,-1/5,4/105,-1/280]; - case 10 - d = [-1/1260,5/504,-5/84,5/21,-5/6,0,5/6,-5/21,5/84,-5/504,1/1260]; - case 12 - d = [1/5544,-1/385,1/56,-5/63,15/56,-6/7,0,6/7,-15/56,5/63,-1/56,1/385,-1/5544]; -end +order = 8; +d = [1/280,-4/105,1/5,-4/5,0,4/5,-1/5,4/105,-1/280]; d = repmat(d,N,1); Q = spdiags(d,-order/2:order/2,N,N);
--- a/+sbp/D1Nonequidistant.m Wed Jul 01 11:15:57 2020 +0200 +++ b/+sbp/D1Nonequidistant.m Wed Jul 01 13:43:32 2020 +0200 @@ -27,50 +27,50 @@ switch option case {'Accurate','accurate','A'} - + [x,h] = sbp.util.accurateBoundaryOptimizedGrid(L,m,order); if order == 4 - [obj.D1,obj.H,obj.x,obj.h] = ... - sbp.implementations.d1_noneq_4(m,L); + [obj.D1,obj.H] = ... + sbp.implementations.d1_noneq_4(m,h); elseif order == 6 - [obj.D1,obj.H,obj.x,obj.h] = ... - sbp.implementations.d1_noneq_6(m,L); + [obj.D1,obj.H] = ... + sbp.implementations.d1_noneq_6(m,h); elseif order == 8 - [obj.D1,obj.H,obj.x,obj.h] = ... - sbp.implementations.d1_noneq_8(m,L); + [obj.D1,obj.H] = ... + sbp.implementations.d1_noneq_8(m,h); elseif order == 10 - [obj.D1,obj.H,obj.x,obj.h] = ... - sbp.implementations.d1_noneq_10(m,L); + [obj.D1,obj.H] = ... + sbp.implementations.d1_noneq_10(m,h); elseif order == 12 - [obj.D1,obj.H,obj.x,obj.h] = ... - sbp.implementations.d1_noneq_12(m,L); + [obj.D1,obj.H] = ... + sbp.implementations.d1_noneq_12(m,h); else error('Invalid operator order %d.',order); end case {'Minimal','minimal','M'} - + [x,h] = sbp.util.minimalBoundaryOptimizedGrid(L,m,order); if order == 4 - [obj.D1,obj.H,obj.x,obj.h] = ... - sbp.implementations.d1_noneq_minimal_4(m,L); + [obj.D1,obj.H] = ... + sbp.implementations.d1_noneq_minimal_4(m,h); elseif order == 6 - [obj.D1,obj.H,obj.x,obj.h] = ... - sbp.implementations.d1_noneq_minimal_6(m,L); + [obj.D1,obj.H] = ... + sbp.implementations.d1_noneq_minimal_6(m,h); elseif order == 8 - [obj.D1,obj.H,obj.x,obj.h] = ... - sbp.implementations.d1_noneq_minimal_8(m,L); + [obj.D1,obj.H] = ... + sbp.implementations.d1_noneq_minimal_8(m,h); elseif order == 10 - [obj.D1,obj.H,obj.x,obj.h] = ... - sbp.implementations.d1_noneq_minimal_10(m,L); + [obj.D1,obj.H] = ... + sbp.implementations.d1_noneq_minimal_10(m,h); elseif order == 12 - [obj.D1,obj.H,obj.x,obj.h] = ... - sbp.implementations.d1_noneq_minimal_12(m,L); + [obj.D1,obj.H] = ... + sbp.implementations.d1_noneq_minimal_12(m,h); else error('Invalid operator order %d.',order); end end - - obj.x = obj.x + x_l; + obj.h = h; + obj.x = x + x_l; obj.e_l = sparse(m,1); obj.e_r = sparse(m,1);