Mercurial > repos > public > sbplib
changeset 1328:3a286d2d1939 feature/D2_boundary_opt
Merge with feature/FMMlabb
| author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
|---|---|
| date | Mon, 14 Feb 2022 11:14:46 +0100 |
| parents | 7ab7d42a5b24 (diff) 74eec7e69b63 (current diff) |
| children | 7df63b17e078 |
| files | |
| diffstat | 22 files changed, 1809 insertions(+), 524 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/boundaryOptimized.m Mon Feb 14 11:14:46 2022 +0100 @@ -0,0 +1,54 @@ +% Creates a Cartesian grid of dimension length(m) +% over the domain xlim, ylim, ... +% The grid is non-equidistant in the boundary regions, +% with node placement based on boundary-optimized SBP operators. +% Examples: +% g = grid.boundaryOptimized([mx, my], xlim, ylim, order, opt) +% g = grid.boundaryOptimized([10, 15], {0,1}, {0,2}, 4) - defaults to 'accurate' stencils +% g = grid.boundaryOptimized([10, 15], {0,1}, {0,2}, 4, 'minimal') +function g = boundaryOptimized(m, varargin) + n = length(m); + + % Check that parameters matches dimensions + matchingParams = false; + if length(varargin) == n+1 % Minimal number of arguments + matchingParams = iscell([varargin{1:n}]) && ... + isfloat([varargin{n+1}]); + elseif length(varargin) == n+2 % Stencil options supplied + matchingParams = iscell([varargin{1:n}]) && ... + isfloat([varargin{n+1}]) && ... + ischar([varargin{n+2}]); + end + assert(matchingParams,'grid:boundaryOptimized:NonMatchingParameters','The number of parameters per dimensions do not match.'); + + % Check that stencil options are passed correctly (if supplied) + if length(varargin) == n+2 % Stencil options supplied + availabe_opts = ["Accurate","accurate","A","Minimal","minimal","M"]; + assert(any(varargin{n+2} == availabe_opts), ... + 'grid:boundaryOptimized:InvalidOption',"The operator option must be 'accurate' or 'minimal.'"); + else %If not passed, populate varargin with default option 'accurate' + varargin(n+2) = {'accurate'}; + end + + % Specify generating function + switch varargin{n+2} + case {'Accurate','accurate','A'} + gridgenerator = @sbp.grid.accurateBoundaryOptimizedGrid; + case {'Minimal','minimal','M'} + gridgenerator = @sbp.grid.minimalBoundaryOptimizedGrid; + end + + X = {}; + h = []; + for i = 1:n + try + [X{i},h(i)] = gridgenerator(varargin{i},m(i),varargin{n+1}); + catch exception % Propagate any errors in the grid generation functions. + msgText = getReport(exception); + error('grid:boundaryOptimized:InvalidParameter',msgText) + end + end + + g = grid.Cartesian(X{:}); + g.h = h; +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/boundaryOptimizedTest.m Mon Feb 14 11:14:46 2022 +0100 @@ -0,0 +1,111 @@ +function tests = boundaryOptimizedTest() + tests = functiontests(localfunctions); +end + +function testErrorInvalidParam(testCase) + in = { + %Invalid order + {[10 10],{0,1},{0,2},3}, + %Invalid grid size + {5, {0,1}, 4}, + {[10 5],{0,1},{0,2},4}, + {[10 5],{0,1},{0,2},6,'M'}, + %Invalid limits + {10,{1},4}, + {[10,10],{0,1},{1},4}, + {[10,10],{1},{1,0},4}, + {10,{1,0},4}, + {[10, 5],{1,0},{0,-1},4}, + }; + + for i = 1:length(in) + testCase.verifyError(@()grid.boundaryOptimized(in{i}{:}),'grid:boundaryOptimized:InvalidParameter',sprintf('in(%d) = %s',i,toString(in{i}))); + end +end + +function testErrorInvalidOption(testCase) + in = { + {[8 8],{0,1},{0,2},4,'acrurate'}, + }; + + for i = 1:length(in) + testCase.verifyError(@()grid.boundaryOptimized(in{i}{:}),'grid:boundaryOptimized:InvalidOption',sprintf('in(%d) = %s',i,toString(in{i}))); + end +end + +function testErrorNonMatchingParam(testCase) + in = { + {[],{1},4}, + {[],{0,1},{0,1},4}, + {[5,5],{0,1},{0,1},{0,1},4}, + {[5,5,4],{0,1},{0,1},4,'accurate'} + {[5,5,4],{0,1},{0,1},{0,1},4,4}, + }; + + for i = 1:length(in) + testCase.verifyError(@()grid.boundaryOptimized(in{i}{:}),'grid:boundaryOptimized:NonMatchingParameters',sprintf('in(%d) = %s',i,toString(in{i}))); + end +end + +% Tests that the expected grid points are obtained for a boundary optimized grid with a 4th order +% accurate stencil and 8th order minimal stencil. +% The boundary grid point distance weights are taken from the D1Nonequidistant operators and +% grid spacing is calculated according to Mattsson et al 2018. +function testCompiles(testCase) + + %% 1D 4th order accurate stencil + % Boundary weights, number of non-equidistantly spaced points for 4th order accurate stencil + bw = [0.0000000000000e+00 6.8764546205559e-01 1.8022115125776e+00]; + n = length(bw)-1; + xi_n = bw(end); + + % Grid points in x-direction. + Lx = 1; + mx = 8; + hx_4 = Lx/(2*xi_n+(mx-2*n-1)); + + bp_l = hx_4*bw; + bp_r = Lx-flip(hx_4*bw); + interior = [hx_4*(xi_n+1) hx_4*(xi_n+2)]; + x_4 = [bp_l interior bp_r]; + + % Boundary weights, number of non-equidistantly spaced points for 8th order minimal stencil + bw = [0.0000000000000e+00, 4.9439570885261e-01, 1.4051531374839e+00]; + n = length(bw)-1; + xi_n = bw(end); + + %% 2D 8th order minimal stencil + % Grid points in x-direction. + hx_8 = Lx/(2*xi_n+(mx-2*n-1)); + + bp_l = hx_8*bw; + bp_r = Lx-flip(hx_8*bw); + interior = [hx_8*(xi_n+1) hx_8*(xi_n+2)]; + x_8 = [bp_l interior bp_r]; + + % Grid points in y-direction. + Ly = 2; + my = 9; + hy = Ly/(2*xi_n+(my-2*n-1)); + + bp_l = hy*bw; + bp_r = Ly-flip(hy*bw); + interior = [hy*(xi_n+1) hy*(xi_n+2) hy*(xi_n+3)]; + y = [bp_l interior bp_r]; + + in = { + {mx, {0,Lx},4}, + {[mx, my],{0,Lx},{0,Ly},8,'M'}, + }; + + out = { + {[x_4'],hx_4} + {[kr(x_8',ones(size(y'))),kr(ones(size(x_8')),y')],[hx_8, hy]} + }; + + for i = 1:length(in) + g = grid.boundaryOptimized(in{i}{:}); + testCase.verifyEqual(g.points(),out{i}{1},'AbsTol', 1e-14, 'RelTol', 1e-14); + testCase.verifyEqual(g.scaling(),out{i}{2},'AbsTol', 1e-14, 'RelTol', 1e-14); + end +end \ No newline at end of file
--- a/+multiblock/+domain/Rectangle.m Fri Mar 13 12:07:00 2020 +0100 +++ b/+multiblock/+domain/Rectangle.m Mon Feb 14 11:14:46 2022 +0100 @@ -93,7 +93,7 @@ N_id = flat_index(m,j,1); S{j} = {S_id,'s'}; N{j} = {N_id,'n'}; - end + end boundaryGroups.E = multiblock.BoundaryGroup(E); boundaryGroups.W = multiblock.BoundaryGroup(W); boundaryGroups.S = multiblock.BoundaryGroup(S); @@ -117,6 +117,20 @@ % Returns a multiblock.Grid given some parameters % ms: cell array of [mx, my] vectors % For same [mx, my] in every block, just input one vector. + % Currently defaults to equidistant grid if varargin is empty. + % If varargin is non-empty, the first argument should supply the grid type, followed by + % additional arguments required to construct the grid. + % Grid types: + % 'equidist' - equidistant grid + % Additional argumets: none + % 'boundaryopt' - boundary optimized grid based on boundary + % optimized SBP operators + % Additional arguments: order, stencil option + % Example: g = getGrid() - the local blocks are 21x21 equidistant grids. + % g = getGrid(ms,) - block i is an equidistant grid with size given by ms{i}. + % g = getGrid(ms,'equidist') - block i is an equidistant grid with size given by ms{i}. + % g = getGrid(ms,'boundaryopt',4,'minimal') - block i is a Cartesian grid with size given by ms{i} + % and nodes placed according to the boundary optimized minimal 4th order SBP operator. function g = getGrid(obj, ms, varargin) default_arg('ms',[21,21]) @@ -129,10 +143,19 @@ ms{i} = m; end end - + if isempty(varargin) || strcmp(varargin{1},'equidist') + gridgenerator = @(m,xlim,ylim)grid.equidistant(m,xlim,ylim); + elseif strcmp(varargin{1},'boundaryopt') + order = varargin{2}; + stenciloption = varargin{3}; + gridgenerator = @(m,xlim,ylim)grid.boundaryOptimized(m,xlim,ylim,... + order,stenciloption); + else + error('No grid type supplied!'); + end grids = cell(1, obj.nBlocks); for i = 1:obj.nBlocks - grids{i} = grid.equidistant(ms{i}, obj.xlims{i}, obj.ylims{i}); + grids{i} = gridgenerator(ms{i}, obj.xlims{i}, obj.ylims{i}); end g = multiblock.Grid(grids, obj.connections, obj.boundaryGroups);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+grid/accurateBoundaryOptimizedGrid.m Mon Feb 14 11:14:46 2022 +0100 @@ -0,0 +1,67 @@ +% Computes the grid points x and grid spacing h used by the boundary optimized SBP operators +% with improved boundary accuracy, presented in +% 'Boundary optimized diagonal-norm SBP operators - Mattsson, Almquist, van der Weide 2018'. +% +% lim - cell array with domain limits +% N - Number of grid points +% order - order of accuracy of sbp operator. +function [x,h] = accurateBoundaryOptimizedGrid(lim,N,order) + assert(iscell(lim) && numel(lim) == 2,'The limit should be cell array with 2 elements.'); + L = lim{2} - lim{1}; + assert(L>0,'Limits must be given in increasing order.'); + %%%% Non-equidistant grid points %%%%% + xb = boundaryPoints(order); + m = length(xb)-1; % Number of non-equidistant points + assert(N-2*(m+1)>=0,'Not enough grid points to contain the boundary region. Requires at least %d points.',2*(m+1)); + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + %%%% 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) ]; + x = x + lim{1}; + %%%%%%%%%%%%%%%%%%%%%%%%% +end +function xb = boundaryPoints(order) + switch order + case 4 + x0 = 0.0000000000000e+00; + x1 = 6.8764546205559e-01; + x2 = 1.8022115125776e+00; + xb = [x0 x1 x2]'; + case 6 + x0 = 0.0000000000000e+00; + x1 = 4.4090263368623e-01; + x2 = 1.2855984345073e+00; + x3 = 2.2638953951239e+00; + xb = [x0 x1 x2 x3]'; + case 8 + x0 = 0.0000000000000e+00; + x1 = 3.8118550247622e-01; + x2 = 1.1899550868338e+00; + x3 = 2.2476300175641e+00; + x4 = 3.3192851303204e+00; + xb = [x0 x1 x2 x3 x4]'; + case 10 + x0 = 0.0000000000000e+00; + x1 = 3.5902433622052e-01; + x2 = 1.1436659188355e+00; + x3 = 2.2144895894456e+00; + x4 = 3.3682742337736e+00; + x5 = 4.4309689056870e+00; + xb = [x0 x1 x2 x3 x4 x5]'; + case 12 + 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; + xb = [x0 x1 x2 x3 x4 x5 x6]'; + otherwise + error('Invalid operator order %d.',order); + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+grid/minimalBoundaryOptimizedGrid.m Mon Feb 14 11:14:46 2022 +0100 @@ -0,0 +1,60 @@ +% Computes the grid points x and grid spacing h used by the boundary optimized SBP operators +% with minimal number of non-equidistant boundary points, presented in +% 'Boundary optimized diagonal-norm SBP operators - Mattsson, Almquist, van der Weide 2018'. +% +% lim - cell array with domain limits +% N - Number of grid points +% order - order of accuracy of sbp operator. +function [x,h] = minimalBoundaryOptimizedGrid(lim,N,order) + assert(iscell(lim) && numel(lim) == 2,'The limit should be a cell array with 2 elements.'); + L = lim{2} - lim{1}; + assert(L>0,'Limits must be given in increasing order.'); + %%%% Non-equidistant grid points %%%%% + xb = boundaryPoints(order); + m = length(xb)-1; % Number of non-equidistant points + assert(N-2*(m+1)>=0,'Not enough grid points to contain the boundary region. Requires at least %d points.',2*(m+1)); + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + %%%% 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) ]; + x = x + lim{1}; + %%%%%%%%%%%%%%%%%%%%%%%%% +end + +function xb = boundaryPoints(order) + switch order + case 4 + x0 = 0.0000000000000e+00; + x1 = 7.7122987842562e-01; + xb = [x0 x1]'; + case 6 + x0 = 0.0000000000000e+00; + x1 = 4.0842950991998e-01; + x2 = 1.1968523189207e+00; + xb = [x0 x1 x2]'; + case 8 + x0 = 0.0000000000000e+00; + x1 = 4.9439570885261e-01; + x2 = 1.4051531374839e+00; + xb = [x0 x1 x2]'; + case 10 + x0 = 0.0000000000000e+00; + x1 = 5.8556160757529e-01; + x2 = 1.7473267488572e+00; + x3 = 3.0000000000000e+00; + xb = [x0 x1 x2 x3]'; + case 12 + x0 = 0.0000000000000e+00; + x1 = 4.6552112904489e-01; + x2 = 1.4647984306493e+00; + x3 = 2.7620429464763e+00; + x4 = 4.0000000000000e+00; + xb = [x0 x1 x2 x3 x4]'; + otherwise + error('Invalid operator order %d.',order); + end +end \ No newline at end of file
--- a/+sbp/+implementations/d1_noneq_10.m Fri Mar 13 12:07:00 2020 +0100 +++ b/+sbp/+implementations/d1_noneq_10.m Mon Feb 14 11:14:46 2022 +0100 @@ -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 Fri Mar 13 12:07:00 2020 +0100 +++ b/+sbp/+implementations/d1_noneq_12.m Mon Feb 14 11:14:46 2022 +0100 @@ -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 Fri Mar 13 12:07:00 2020 +0100 +++ b/+sbp/+implementations/d1_noneq_4.m Mon Feb 14 11:14:46 2022 +0100 @@ -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 Fri Mar 13 12:07:00 2020 +0100 +++ b/+sbp/+implementations/d1_noneq_6.m Mon Feb 14 11:14:46 2022 +0100 @@ -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 Fri Mar 13 12:07:00 2020 +0100 +++ b/+sbp/+implementations/d1_noneq_8.m Mon Feb 14 11:14:46 2022 +0100 @@ -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 Fri Mar 13 12:07:00 2020 +0100 +++ b/+sbp/+implementations/d1_noneq_minimal_10.m Mon Feb 14 11:14:46 2022 +0100 @@ -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 Fri Mar 13 12:07:00 2020 +0100 +++ b/+sbp/+implementations/d1_noneq_minimal_12.m Mon Feb 14 11:14:46 2022 +0100 @@ -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 Fri Mar 13 12:07:00 2020 +0100 +++ b/+sbp/+implementations/d1_noneq_minimal_4.m Mon Feb 14 11:14:46 2022 +0100 @@ -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 Fri Mar 13 12:07:00 2020 +0100 +++ b/+sbp/+implementations/d1_noneq_minimal_6.m Mon Feb 14 11:14:46 2022 +0100 @@ -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 Fri Mar 13 12:07:00 2020 +0100 +++ b/+sbp/+implementations/d1_noneq_minimal_8.m Mon Feb 14 11:14:46 2022 +0100 @@ -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);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/d2_noneq_variable_10.m Mon Feb 14 11:14:46 2022 +0100 @@ -0,0 +1,320 @@ +function [H, HI, D1, D2, DI] = d2_noneq_variable_10(N, h, options) + % N: Number of grid points + % h: grid spacing + % options: struct containing options for constructing the operator + % current options are: + % options.stencil_type ('minimal','nonminimal','wide') + % options.AD ('upwind', 'op') + + % BP: Number of boundary points + % order: Accuracy of interior stencil + BP = 10; + order = 10; + + %%%% Norm matrix %%%%%%%% + P = zeros(BP, 1); + P0 = 1.0000000000000e-01; + P1 = 5.8980851260667e-01; + P2 = 9.5666820955973e-01; + P3 = 1.1500297411596e+00; + P4 = 1.1232986993248e+00; + P5 = 1.0123020150951e+00; + P6 = 9.9877122702527e-01; + P7 = 1.0000873322761e+00; + P8 = 1.0000045540888e+00; + P9 = 9.9999861455083e-01; + + for i = 0:BP - 1 + P(i + 1) = eval(['P' num2str(i)]); + end + + Hv = ones(N, 1); + Hv(1:BP) = P; + Hv(end - BP + 1:end) = flip(P); + Hv = h * Hv; + H = spdiags(Hv, 0, N, N); + HI = spdiags(1 ./ Hv, 0, N, N); + %%%%%%%%%%%%%%%%%%%%%%%%% + + %%%% Q matrix %%%%%%%%%%% + + % interior stencil + 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); + + % Boundaries + Q0_0 = -5.0000000000000e-01; + Q0_1 = 6.7548747038002e-01; + Q0_2 = -2.6691978151546e-01; + Q0_3 = 1.4438714982130e-01; + Q0_4 = -7.7273673750760e-02; + Q0_5 = 2.5570078343005e-02; + Q0_6 = 4.2808774693299e-03; + Q0_7 = -8.2902108933389e-03; + Q0_8 = 3.2031176427908e-03; + Q0_9 = -4.4502749689556e-04; + Q0_10 = 0.0000000000000e+00; + Q0_11 = 0.0000000000000e+00; + Q0_12 = 0.0000000000000e+00; + Q0_13 = 0.0000000000000e+00; + Q0_14 = 0.0000000000000e+00; + Q1_0 = -6.7548747038002e-01; + Q1_1 = 0.0000000000000e+00; + Q1_2 = 9.5146052715180e-01; + Q1_3 = -4.2442349882626e-01; + Q1_4 = 2.1538865145190e-01; + Q1_5 = -7.1939778160350e-02; + Q1_6 = -8.2539187832840e-03; + Q1_7 = 1.9930661669090e-02; + Q1_8 = -7.7433256989613e-03; + Q1_9 = 1.0681515760869e-03; + Q1_10 = 0.0000000000000e+00; + Q1_11 = 0.0000000000000e+00; + Q1_12 = 0.0000000000000e+00; + Q1_13 = 0.0000000000000e+00; + Q1_14 = 0.0000000000000e+00; + Q2_0 = 2.6691978151546e-01; + Q2_1 = -9.5146052715180e-01; + Q2_2 = 0.0000000000000e+00; + Q2_3 = 9.6073770842387e-01; + Q2_4 = -3.9378595264609e-01; + Q2_5 = 1.3302097358959e-01; + Q2_6 = 8.1200458151489e-05; + Q2_7 = -2.3849770528789e-02; + Q2_8 = 9.6600442856829e-03; + Q2_9 = -1.3234579460680e-03; + Q2_10 = 0.0000000000000e+00; + Q2_11 = 0.0000000000000e+00; + Q2_12 = 0.0000000000000e+00; + Q2_13 = 0.0000000000000e+00; + Q2_14 = 0.0000000000000e+00; + Q3_0 = -1.4438714982130e-01; + Q3_1 = 4.2442349882626e-01; + Q3_2 = -9.6073770842387e-01; + Q3_3 = 0.0000000000000e+00; + Q3_4 = 9.1551097634196e-01; + Q3_5 = -2.8541713079648e-01; + Q3_6 = 4.1398809121293e-02; + Q3_7 = 1.7256059167927e-02; + Q3_8 = -9.4349194803610e-03; + Q3_9 = 1.3875650645663e-03; + Q3_10 = 0.0000000000000e+00; + Q3_11 = 0.0000000000000e+00; + Q3_12 = 0.0000000000000e+00; + Q3_13 = 0.0000000000000e+00; + Q3_14 = 0.0000000000000e+00; + Q4_0 = 7.7273673750760e-02; + Q4_1 = -2.1538865145190e-01; + Q4_2 = 3.9378595264609e-01; + Q4_3 = -9.1551097634196e-01; + Q4_4 = 0.0000000000000e+00; + Q4_5 = 8.3519401865051e-01; + Q4_6 = -2.0586492924974e-01; + Q4_7 = 3.1230261235901e-02; + Q4_8 = -2.0969453466651e-04; + Q4_9 = -5.0965470499782e-04; + Q4_10 = 0.0000000000000e+00; + Q4_11 = 0.0000000000000e+00; + Q4_12 = 0.0000000000000e+00; + Q4_13 = 0.0000000000000e+00; + Q4_14 = 0.0000000000000e+00; + Q5_0 = -2.5570078343005e-02; + Q5_1 = 7.1939778160350e-02; + Q5_2 = -1.3302097358959e-01; + Q5_3 = 2.8541713079648e-01; + Q5_4 = -8.3519401865051e-01; + Q5_5 = 0.0000000000000e+00; + Q5_6 = 8.1046389580138e-01; + Q5_7 = -2.1879194972141e-01; + Q5_8 = 5.2977237804899e-02; + Q5_9 = -9.0146730522360e-03; + Q5_10 = 7.9365079365079e-04; + Q5_11 = 0.0000000000000e+00; + Q5_12 = 0.0000000000000e+00; + Q5_13 = 0.0000000000000e+00; + Q5_14 = 0.0000000000000e+00; + Q6_0 = -4.2808774693299e-03; + Q6_1 = 8.2539187832840e-03; + Q6_2 = -8.1200458151489e-05; + Q6_3 = -4.1398809121293e-02; + Q6_4 = 2.0586492924974e-01; + Q6_5 = -8.1046389580138e-01; + Q6_6 = 0.0000000000000e+00; + Q6_7 = 8.2787884456005e-01; + Q6_8 = -2.3582460382545e-01; + Q6_9 = 5.9178678209520e-02; + Q6_10 = -9.9206349206349e-03; + Q6_11 = 7.9365079365079e-04; + Q6_12 = 0.0000000000000e+00; + Q6_13 = 0.0000000000000e+00; + Q6_14 = 0.0000000000000e+00; + Q7_0 = 8.2902108933389e-03; + Q7_1 = -1.9930661669090e-02; + Q7_2 = 2.3849770528789e-02; + Q7_3 = -1.7256059167927e-02; + Q7_4 = -3.1230261235901e-02; + Q7_5 = 2.1879194972141e-01; + Q7_6 = -8.2787884456005e-01; + Q7_7 = 0.0000000000000e+00; + Q7_8 = 8.3301028859275e-01; + Q7_9 = -2.3804321850015e-01; + Q7_10 = 5.9523809523809e-02; + Q7_11 = -9.9206349206349e-03; + Q7_12 = 7.9365079365079e-04; + Q7_13 = 0.0000000000000e+00; + Q7_14 = 0.0000000000000e+00; + Q8_0 = -3.2031176427908e-03; + Q8_1 = 7.7433256989613e-03; + Q8_2 = -9.6600442856829e-03; + Q8_3 = 9.4349194803610e-03; + Q8_4 = 2.0969453466651e-04; + Q8_5 = -5.2977237804899e-02; + Q8_6 = 2.3582460382545e-01; + Q8_7 = -8.3301028859275e-01; + Q8_8 = 0.0000000000000e+00; + Q8_9 = 8.3333655748509e-01; + Q8_10 = -2.3809523809524e-01; + Q8_11 = 5.9523809523809e-02; + Q8_12 = -9.9206349206349e-03; + Q8_13 = 7.9365079365079e-04; + Q8_14 = 0.0000000000000e+00; + Q9_0 = 4.4502749689556e-04; + Q9_1 = -1.0681515760869e-03; + Q9_2 = 1.3234579460680e-03; + Q9_3 = -1.3875650645663e-03; + Q9_4 = 5.0965470499782e-04; + Q9_5 = 9.0146730522360e-03; + Q9_6 = -5.9178678209520e-02; + Q9_7 = 2.3804321850015e-01; + Q9_8 = -8.3333655748509e-01; + Q9_9 = 0.0000000000000e+00; + Q9_10 = 8.3333333333333e-01; + Q9_11 = -2.3809523809524e-01; + Q9_12 = 5.9523809523809e-02; + Q9_13 = -9.9206349206349e-03; + Q9_14 = 7.9365079365079e-04; + + for i = 1:BP + + for j = 1:BP + Q(i, j) = eval(['Q' num2str(i - 1) '_' num2str(j - 1)]); + Q(N + 1 - i, N + 1 - j) = -eval(['Q' num2str(i - 1) '_' num2str(j - 1)]); + end + + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%% + + %%% Undivided difference operators %%%% + % Closed with zeros at the first boundary nodes. + m = N; + + DD_5 = (-diag(ones(m - 3, 1), -3) + 5 * diag(ones(m - 2, 1), -2) - 10 * diag(ones(m - 1, 1), -1) + 10 * diag(ones(m, 1), 0) - 5 * diag(ones(m - 1, 1), 1) + diag(ones(m - 2, 1), 2)); + DD_5(1:8, 1:10) = [0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0; -0.88425676349344827581e1 0.18735837038089672208e2 -0.17076524806228533350e2 0.10664693462778093925e2 -0.43403898038449865885e1 0.85895174414023656383e0 0 0 0 0; 0 -0.13262411219221127021e1 0.45552739083119183335e1 -0.73424543813164879330e1 0.70876331528332015948e1 -0.38059884697709954324e1 0.83177691186447613913e0 0 0 0; 0 0 -0.67600766166626540709e0 0.32310531977344450754e1 -0.69639522132755183916e1 0.77488870404680024829e1 -0.42187263911386203939e1 0.87874602787795663432e0 0 0; 0 0 0 -0.66326118352702166971e0 0.38132489024062224761e1 -0.84909798376351149047e1 0.90434791287189283383e1 -0.46461964978830115989e1 0.94370948791999735890e0 0; 0 0 0 0 -0.86903681018048694197e0 0.47050202961852482685e1 -0.96960545214617434798e1 0.97952957162584740145e1 -0.49228410242754295559e1 0.98761634347393769455e0; ]; + DD_5(m - 6:m, m - 9:m) = [-0.98761634347393769455e0 0.49228410242754295559e1 -0.97952957162584740145e1 0.96960545214617434798e1 -0.47050202961852482685e1 0.86903681018048694197e0 0 0 0 0; 0 -0.94370948791999735890e0 0.46461964978830115989e1 -0.90434791287189283383e1 0.84909798376351149047e1 -0.38132489024062224761e1 0.66326118352702166971e0 0 0 0; 0 0 -0.87874602787795663432e0 0.42187263911386203939e1 -0.77488870404680024829e1 0.69639522132755183916e1 -0.32310531977344450754e1 0.67600766166626540709e0 0 0; 0 0 0 -0.83177691186447613913e0 0.38059884697709954324e1 -0.70876331528332015948e1 0.73424543813164879330e1 -0.45552739083119183335e1 0.13262411219221127021e1 0; 0 0 0 0 -0.85895174414023656383e0 0.43403898038449865885e1 -0.10664693462778093925e2 0.17076524806228533350e2 -0.18735837038089672208e2 0.88425676349344827581e1; 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0; ]; + DD_5 = sparse(DD_5); + + DD_6 = (diag(ones(m - 3, 1), 3) - 6 * diag(ones(m - 2, 1), 2) + 15 * diag(ones(m - 1, 1), 1) - 20 * diag(ones(m, 1), 0) + 15 * diag(ones(m - 1, 1), -1) - 6 * diag(ones(m - 2, 1), -2) + diag(ones(m - 3, 1), -3)); + DD_6(1:8, 1:11) = [0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0; 0.97690498198305661286e1 -0.22164087301995739032e2 0.23898275711233304083e2 -0.19893851160044639280e2 0.12625396854743289626e2 -0.51537104648414193830e1 0.91892654107463785790e0 0 0 0 0; 0 0.13105269062480722931e1 -0.51692977530964314903e1 0.10448225399376470222e2 -0.13885092532071720368e2 0.11417965409312986297e2 -0.49906614711868568348e1 0.86833404141747988109e0 0 0 0; 0 0 0.64511698871199831925e0 -0.37163607887886698482e1 0.10284728893981919286e2 -0.15497774080936004966e2 0.12656179173415861182e2 -0.52724761672677398059e1 0.90058598088263583338e0 0 0; 0 0 0 0.64016413450696574250e0 -0.45192323253005276829e1 0.12736469756452672357e2 -0.18086958257437856677e2 0.13938589493649034797e2 -0.56622569275199841534e1 0.95322412564969561677e0 0; 0 0 0 0 0.86005005088560527649e0 -0.56460243554222979222e1 0.14544081782192615220e2 -0.19590591432516948029e2 0.14768523072826288668e2 -0.59256980608436261673e1 0.98965894287836295484e0; ]; + DD_6(m - 7:m, m - 10:m) = [0.98965894287836295484e0 -0.59256980608436261673e1 0.14768523072826288668e2 -0.19590591432516948029e2 0.14544081782192615220e2 -0.56460243554222979222e1 0.86005005088560527649e0 0 0 0 0; 0 0.95322412564969561677e0 -0.56622569275199841534e1 0.13938589493649034797e2 -0.18086958257437856677e2 0.12736469756452672357e2 -0.45192323253005276829e1 0.64016413450696574250e0 0 0 0; 0 0 0.90058598088263583338e0 -0.52724761672677398059e1 0.12656179173415861182e2 -0.15497774080936004966e2 0.10284728893981919286e2 -0.37163607887886698482e1 0.64511698871199831925e0 0 0; 0 0 0 0.86833404141747988109e0 -0.49906614711868568348e1 0.11417965409312986297e2 -0.13885092532071720368e2 0.10448225399376470222e2 -0.51692977530964314903e1 0.13105269062480722931e1 0; 0 0 0 0 0.91892654107463785790e0 -0.51537104648414193830e1 0.12625396854743289626e2 -0.19893851160044639280e2 0.23898275711233304083e2 -0.22164087301995739032e2 0.97690498198305661286e1; 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0; ]; + DD_6 = sparse(DD_6); + + DD_7 = (-diag(ones(m - 4, 1), -4) + 7 * diag(ones(m - 3, 1), -3) - 21 * diag(ones(m - 2, 1), -2) + 35 * diag(ones(m - 1, 1), -1) - 35 * diag(ones(m, 1), 0) + 21 * diag(ones(m - 1, 1), 1) - 7 * diag(ones(m - 2, 1), 2) + diag(ones(m - 3, 1), 3)); + DD_7(1:9, 1:12) = [0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0; -0.10633444157744515342e2 0.25551717302255079868e2 -0.31639558087487298930e2 0.33026832975062979647e2 -0.28856215669710369771e2 0.18037986626944967840e2 -0.64324857875224650053e1 0.94516679820162169311e0 0 0 0 0; 0 -0.12971946051930944772e1 0.57552633215463132945e1 -0.14020486493241325087e2 0.23923936099960072962e2 -0.26641919288396968027e2 0.17467315149153998922e2 -0.60783382899223591677e1 0.89142410609336157976e0 0 0 0; 0 0 -0.61968315701047317846e0 0.41847682905586435755e1 -0.14220312881453249947e2 0.27121104641638008690e2 -0.29531084737970342758e2 0.18453666585437089321e2 -0.63041018661784508337e1 0.91564312497877513070e0 0 0; 0 0 0 -0.62096054671195295664e0 0.52179151332917540372e1 -0.17831057659033741300e2 0.31652176950516249184e2 -0.32523375485181081192e2 0.19817899246319944537e2 -0.66725688795478693174e1 0.95997124034669700791e0 0; 0 0 0 0 -0.85241549236735026150e0 0.65870284146593475759e1 -0.20361714495069661308e2 0.34283535006904659051e2 -0.34459887169928006891e2 0.20739943212952691585e2 -0.69276126001485406839e1 0.99112312299686093167e0; ]; + DD_7(m - 7:m, m - 11:m) = [-0.99112312299686093167e0 0.69276126001485406839e1 -0.20739943212952691585e2 0.34459887169928006891e2 -0.34283535006904659051e2 0.20361714495069661308e2 -0.65870284146593475759e1 0.85241549236735026150e0 0 0 0 0; 0 -0.95997124034669700791e0 0.66725688795478693174e1 -0.19817899246319944537e2 0.32523375485181081192e2 -0.31652176950516249184e2 0.17831057659033741300e2 -0.52179151332917540372e1 0.62096054671195295664e0 0 0 0; 0 0 -0.91564312497877513070e0 0.63041018661784508337e1 -0.18453666585437089321e2 0.29531084737970342758e2 -0.27121104641638008690e2 0.14220312881453249947e2 -0.41847682905586435755e1 0.61968315701047317846e0 0 0; 0 0 0 -0.89142410609336157976e0 0.60783382899223591677e1 -0.17467315149153998922e2 0.26641919288396968027e2 -0.23923936099960072962e2 0.14020486493241325087e2 -0.57552633215463132945e1 0.12971946051930944772e1 0; 0 0 0 0 -0.94516679820162169311e0 0.64324857875224650053e1 -0.18037986626944967840e2 0.28856215669710369771e2 -0.33026832975062979647e2 0.31639558087487298930e2 -0.25551717302255079868e2 0.10633444157744515342e2; 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0; ]; + DD_7 = sparse(DD_7); + + DD_8 = (diag(ones(m - 4, 1), 4) - 8 * diag(ones(m - 3, 1), 3) + 28 * diag(ones(m - 2, 1), 2) - 56 * diag(ones(m - 1, 1), 1) + 70 * diag(ones(m, 1), 0) - 56 * diag(ones(m - 1, 1), -1) + 28 * diag(ones(m - 2, 1), -2) - 8 * diag(ones(m - 3, 1), -3) + diag(ones(m - 4, 1), -4)); + DD_8(1:9, 1:13) = [0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0; 0.11447706798618566234e2 -0.28904884139025747734e2 0.40258353260409964223e2 -0.50649997399180126990e2 0.56821824921675878869e2 -0.48101297671853247575e2 0.25729943150089860021e2 -0.75613343856129735449e1 0.95968546487782649665e0 0 0 0 0; 0 0.12856328177474934636e1 -0.63181271116961104382e1 0.18042992864608611563e2 -0.37804272468074727719e2 0.53283838576793936053e2 -0.46579507064410663791e2 0.24313353159689436671e2 -0.71313928487468926381e1 0.90748207408891683663e0 0 0 0; 0 0 0.59820007352804876260e0 -0.46391245450012267365e1 0.18764346418211454409e2 -0.43393767426620813904e2 0.59062169475940685515e2 -0.49209777561165571522e2 0.25216407464713803335e2 -0.73251449998302010456e1 0.92669110022382118702e0 0 0; 0 0 0 0.60460011916248255841e0 -0.59103958199322344450e1 0.23774743545378321733e2 -0.50643483120825998695e2 0.65046750970362162385e2 -0.52847731323519852099e2 0.26690275518191477270e2 -0.76797699227735760633e1 0.96501003395721735560e0 0; 0 0 0 0 0.84578719850252712660e0 -0.75280324738963972296e1 0.27148952660092881743e2 -0.54853656011047454481e2 0.68919774339856013782e2 -0.55306515234540510895e2 0.27710450400594162736e2 -0.79289849839748874534e1 0.99222410441366467138e0; ]; + DD_8(m - 8:m, m - 12:m) = [0.99222410441366467138e0 -0.79289849839748874534e1 0.27710450400594162736e2 -0.55306515234540510895e2 0.68919774339856013782e2 -0.54853656011047454481e2 0.27148952660092881743e2 -0.75280324738963972296e1 0.84578719850252712660e0 0 0 0 0; 0 0.96501003395721735560e0 -0.76797699227735760633e1 0.26690275518191477270e2 -0.52847731323519852099e2 0.65046750970362162385e2 -0.50643483120825998695e2 0.23774743545378321733e2 -0.59103958199322344450e1 0.60460011916248255841e0 0 0 0; 0 0 0.92669110022382118702e0 -0.73251449998302010456e1 0.25216407464713803335e2 -0.49209777561165571522e2 0.59062169475940685515e2 -0.43393767426620813904e2 0.18764346418211454409e2 -0.46391245450012267365e1 0.59820007352804876260e0 0 0; 0 0 0 0.90748207408891683663e0 -0.71313928487468926381e1 0.24313353159689436671e2 -0.46579507064410663791e2 0.53283838576793936053e2 -0.37804272468074727719e2 0.18042992864608611563e2 -0.63181271116961104382e1 0.12856328177474934636e1 0; 0 0 0 0 0.95968546487782649665e0 -0.75613343856129735449e1 0.25729943150089860021e2 -0.48101297671853247575e2 0.56821824921675878869e2 -0.50649997399180126990e2 0.40258353260409964223e2 -0.28904884139025747734e2 0.11447706798618566234e2; 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0; ]; + DD_8 = sparse(DD_8); + + DD_9 = (-diag(ones(m - 5, 1), -5) + 9 * diag(ones(m - 4, 1), -4) - 36 * diag(ones(m - 3, 1), -3) + 84 * diag(ones(m - 2, 1), -2) - 126 * diag(ones(m - 1, 1), -1) + 126 * diag(ones(m, 1), 0) - 84 * diag(ones(m - 1, 1), 1) + 36 * diag(ones(m - 2, 1), 2) - 9 * diag(ones(m - 3, 1), 3) + diag(ones(m - 4, 1), 4)); + DD_9(1:10, 1:14) = [0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0; -0.12220346479758703402e2 0.32228164479135716142e2 -0.49720065159556828467e2 0.73329283892516285142e2 -0.10101269332559123572e3 0.10822791976166980704e3 -0.77189829450269580063e2 0.34026004735258380952e2 -0.86371691839004384699e1 0.96873073049659684460e0 0 0 0 0; 0 -0.12754371756934054072e1 0.68614776237194576426e1 -0.22502238094968228382e2 0.56120004490560654424e2 -0.95910909438229084895e2 0.10480389089492399353e3 -0.72940059479068310012e2 0.32091267819361016871e2 -0.81673386668002515297e1 0.91934202619415775756e0 0 0 0; 0 0 -0.57969473693003825890e0 0.50815098898235167347e1 -0.23911428372444897689e2 0.65090651139931220857e2 -0.10631190505669323393e3 0.11072199951262253592e3 -0.75649222394141410004e2 0.32963152499235904705e2 -0.83402199020143906832e1 0.93515742061079234147e0 0 0; 0 0 0 -0.59039909771982400264e0 0.65974918490578446834e1 -0.30567527415486413657e2 0.75965224681238998042e2 -0.11708415174665189229e3 0.11890739547791966722e3 -0.80070826554574431809e2 0.34558964652481092285e2 -0.86850903056149562004e1 0.96891845934991572940e0 0; 0 0 0 0 -0.83993614063969059233e0 0.84690365331334468834e1 -0.34905796277262276527e2 0.82280484016571181722e2 -0.12405559381174082481e3 0.12443965927771614951e3 -0.83131351201782488207e2 0.35680432427886993540e2 -0.89300169397229820424e1 0.99308211584049051802e0; ]; + DD_9(m - 8:m, m - 13:m) = [-0.99308211584049051802e0 0.89300169397229820424e1 -0.35680432427886993540e2 0.83131351201782488207e2 -0.12443965927771614951e3 0.12405559381174082481e3 -0.82280484016571181722e2 0.34905796277262276527e2 -0.84690365331334468834e1 0.83993614063969059233e0 0 0 0 0; 0 -0.96891845934991572940e0 0.86850903056149562004e1 -0.34558964652481092285e2 0.80070826554574431809e2 -0.11890739547791966722e3 0.11708415174665189229e3 -0.75965224681238998042e2 0.30567527415486413657e2 -0.65974918490578446834e1 0.59039909771982400264e0 0 0 0; 0 0 -0.93515742061079234147e0 0.83402199020143906832e1 -0.32963152499235904705e2 0.75649222394141410004e2 -0.11072199951262253592e3 0.10631190505669323393e3 -0.65090651139931220857e2 0.23911428372444897689e2 -0.50815098898235167347e1 0.57969473693003825890e0 0 0; 0 0 0 -0.91934202619415775756e0 0.81673386668002515297e1 -0.32091267819361016871e2 0.72940059479068310012e2 -0.10480389089492399353e3 0.95910909438229084895e2 -0.56120004490560654424e2 0.22502238094968228382e2 -0.68614776237194576426e1 0.12754371756934054072e1 0; 0 0 0 0 -0.96873073049659684460e0 0.86371691839004384699e1 -0.34026004735258380952e2 0.77189829450269580063e2 -0.10822791976166980704e3 0.10101269332559123572e3 -0.73329283892516285142e2 0.49720065159556828467e2 -0.32228164479135716142e2 0.12220346479758703402e2; 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0; ]; + DD_9 = sparse(DD_9); + + DD_10 = (diag(ones(m - 5, 1), 5) - 10 * diag(ones(m - 4, 1), 4) + 45 * diag(ones(m - 3, 1), 3) - 120 * diag(ones(m - 2, 1), 2) + 210 * diag(ones(m - 1, 1), 1) - 252 * diag(ones(m, 1), 0) + 210 * diag(ones(m - 1, 1), -1) - 120 * diag(ones(m - 2, 1), -2) + 45 * diag(ones(m - 3, 1), -3) - 10 * diag(ones(m - 4, 1), -4) + diag(ones(m - 5, 1), -5)); + DD_10(1:10, 1:15) = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0.12957678688124699986e2 -0.35525089722887367966e2 0.59995471673283665865e2 -0.10161365491269611714e3 0.16661352548983481943e3 -0.21645583952333961409e3 0.19297457362567395016e3 -0.11342001578419460317e3 0.43185845919502192349e2 -0.96873073049659684460e1 0.97481185166434302666e0 0 0 0 0; 0 0.12663266431786644877e1 -0.73880195719183439765e1 0.27386715439651102593e2 -0.79459763018974764721e2 0.15985151573038180816e3 -0.20960778178984798706e3 0.18235014869767077503e3 -0.10697089273120338957e3 0.40836693334001257648e2 -0.91934202619415775756e1 0.92847752900245498778e0 0 0 0; 0 0 0.56350506801536116666e0 -0.55135043604652958562e1 0.29656869502625787574e2 -0.92986644485616029795e2 0.17718650842782205655e3 -0.22144399902524507185e3 0.18912305598535352501e3 -0.10987717499745301568e3 0.41701099510071953416e2 -0.93515742061079234147e1 0.94185858099865288757e0 0 0; 0 0 0 0.57788899624281661128e0 -0.72798346274475049971e1 0.38209409269358017071e2 -0.10852174954462714006e3 0.19514025291108648715e3 -0.23781479095583933444e3 0.20017706638643607952e3 -0.11519654884160364095e3 0.43425451528074781002e2 -0.96891845934991572940e1 0.97203947181859638306e0 0; 0 0 0 0 0.83470299758184640199e0 -0.94100405923704965371e1 0.43632245346577845659e2 -0.11754354859510168817e3 0.20675932301956804135e3 -0.24887931855543229903e3 0.20782837800445622052e3 -0.11893477475962331180e3 0.44650084698614910212e2 -0.99308211584049051802e1 0.99376959413383658153e0; ]; + DD_10(m - 9:m, m - 14:m) = [0.99376959413383658153e0 -0.99308211584049051802e1 0.44650084698614910212e2 -0.11893477475962331180e3 0.20782837800445622052e3 -0.24887931855543229903e3 0.20675932301956804135e3 -0.11754354859510168817e3 0.43632245346577845659e2 -0.94100405923704965371e1 0.83470299758184640199e0 0 0 0 0; 0 0.97203947181859638306e0 -0.96891845934991572940e1 0.43425451528074781002e2 -0.11519654884160364095e3 0.20017706638643607952e3 -0.23781479095583933444e3 0.19514025291108648715e3 -0.10852174954462714006e3 0.38209409269358017071e2 -0.72798346274475049971e1 0.57788899624281661128e0 0 0 0; 0 0 0.94185858099865288757e0 -0.93515742061079234147e1 0.41701099510071953416e2 -0.10987717499745301568e3 0.18912305598535352501e3 -0.22144399902524507185e3 0.17718650842782205655e3 -0.92986644485616029795e2 0.29656869502625787574e2 -0.55135043604652958562e1 0.56350506801536116666e0 0 0; 0 0 0 0.92847752900245498778e0 -0.91934202619415775756e1 0.40836693334001257648e2 -0.10697089273120338957e3 0.18235014869767077503e3 -0.20960778178984798706e3 0.15985151573038180816e3 -0.79459763018974764721e2 0.27386715439651102593e2 -0.73880195719183439765e1 0.12663266431786644877e1 0; 0 0 0 0 0.97481185166434302666e0 -0.96873073049659684460e1 0.43185845919502192349e2 -0.11342001578419460317e3 0.19297457362567395016e3 -0.21645583952333961409e3 0.16661352548983481943e3 -0.10161365491269611714e3 0.59995471673283665865e2 -0.35525089722887367966e2 0.12957678688124699986e2; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; ]; + DD_10 = sparse(DD_10); + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + %%%% Difference operator %% + D1 = H \ Q; + + % Helper functions for constructing D2(c) + % TODO: Consider changing sparse(diag(...)) to spdiags(....) + + % Minimal 11 point stencil width + function D2 = D2_fun_minimal(c) + % Here we add variable diffusion + C1 = sparse(diag(c)); + C2 = 1/2 * diag(ones(m - 1, 1), -1) + 1/2 * diag(ones(m, 1), 0); C2(1, 2) = 1/2; + C3 = 2/7 * diag(ones(m - 1, 1), -1) + 2/7 * diag(ones(m - 1, 1), 1) + 3/7 * diag(ones(m, 1), 0); C3(1, 3) = 2/7; C3(m, m - 2) = 2/7; + C4 = 1/5 * diag(ones(m - 2, 1), -2) + 3/10 * diag(ones(m - 1, 1), -1) + 1/5 * diag(ones(m - 1, 1), 1) + 3/10 * diag(ones(m, 1), 0); C4(2, 4) = 1/5; C4(1, 3) = 1/5; C4(1, 4) = 3/10; C4(m, m - 2) = 1/5; C4(m, m - 3) = 1/5; + C5 = 1/5 * diag(ones(m - 2, 1), -2) + 1/5 * diag(ones(m - 2, 1), 2) + 1/5 * diag(ones(m - 1, 1), -1) + 1/5 * diag(ones(m - 1, 1), 1) + 1/5 * diag(ones(m, 1), 0); C5(2, 4) = 1/5; C5(1, 4) = 1/5; C5(1, 5) = 1/5; C5(2, 5) = 1/5; C5(m - 1, m - 4) = 1/5; C5(m, m - 4) = 1/5; C5(m, m - 3) = 1/5; + + C2 = sparse(diag(C2 * c)); + C3 = sparse(diag(C3 * c)); + C4 = sparse(diag(C4 * c)); + C5 = sparse(diag(C5 * c)); + + % Remainder term added to wide second derivative operator + R = (1/1587600 / h) * transpose(DD_10) * C1 * DD_10 + (1/317520 / h) * transpose(DD_9) * C2 * DD_9 + (1/60480 / h) * transpose(DD_8) * C3 * DD_8 + (1/10584 / h) * transpose(DD_7) * C4 * DD_7 + (1/1512 / h) * transpose(DD_6) * C5 * DD_6; + D2 = D1 * C1 * D1 - H \ R; + end + + % Few additional grid point in interior stencil cmp. to minimal + function D2 = D2_fun_nonminimal(c) + C1 = sparse(diag(c)); + C2 = 1/2 * diag(ones(m - 1, 1), -1) + 1/2 * diag(ones(m, 1), 0); C2(1, 2) = 1/2; + C3 = 2/7 * diag(ones(m - 1, 1), -1) + 2/7 * diag(ones(m - 1, 1), 1) + 3/7 * diag(ones(m, 1), 0); C3(1, 3) = 2/7; C3(m, m - 2) = 2/7; + C4 = 1/5 * diag(ones(m - 2, 1), -2) + 3/10 * diag(ones(m - 1, 1), -1) + 1/5 * diag(ones(m - 1, 1), 1) + 3/10 * diag(ones(m, 1), 0); C4(2, 4) = 1/5; C4(1, 3) = 1/5; C4(1, 4) = 3/10; C4(m, m - 2) = 1/5; C4(m, m - 3) = 1/5; + + C2 = sparse(diag(C2 * c)); + C3 = sparse(diag(C3 * c)); + C4 = sparse(diag(C4 * c)); + + % Remainder term added to wide second derivative operator + R = (1/1587600 / h) * transpose(DD_10) * C1 * DD_10 + (1/317520 / h) * transpose(DD_9) * C2 * DD_9 + (1/60480 / h) * transpose(DD_8) * C3 * DD_8 + (1/10584 / h) * transpose(DD_7) * C4 * DD_7; + D2 = D1 * C1 * D1 - H \ R; + end + + % Wide stencil + function D2 = D2_fun_wide(c) + % Here we add variable diffusion + C1 = sparse(diag(c)); + D2 = D1 * C1 * D1; + end + + switch options.stencil_width + case 'minimal' + D2 = @D2_fun_minimal; + case 'nonminimal' + D2 = @D2_fun_nonminimal; + case 'wide' + D2 = @D2_fun_wide; + otherwise + error('No option %s for stencil width', options.stencil_width) + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%% + + %%%% Artificial dissipation operator %%% + switch options.AD + case 'upwind' + % This is the choice that yield 9th order Upwind + DI = H \ (transpose(DD_5) * DD_5) * (-1/1260); + case 'op' + % This choice will preserve the order of the underlying + % Non-dissipative D1 SBP operator + DI = H \ (transpose(DD_6) * DD_6) * (-1 / (5 * 1260)); + % Notice that you can use any negative number instead of (-1/(5*1260)) + otherwise + error("Artificial dissipation options '%s' not implemented.", option.AD) + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%% +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/d2_noneq_variable_12.m Mon Feb 14 11:14:46 2022 +0100 @@ -0,0 +1,401 @@ +function [H, HI, D1, D2, DI] = d2_noneq_variable_12(N, h, options) + % N: Number of grid points + % h: grid spacing + % options: struct containing options for constructing the operator + % current options are: + % options.stencil_type ('minimal','nonminimal','wide') + % options.AD ('upwind', 'op') + + % BP: Number of boundary points + % order: Accuracy of interior stencil + BP = 12; + order = 12; + + %%%% Norm matrix %%%%%%%% + P = zeros(BP, 1); + P0 = 1.0000000000011e-01; + P1 = 5.9616216757547e-01; + P2 = 9.9065699844442e-01; + P3 = 1.2512548713913e+00; + P4 = 1.3316678989403e+00; + P5 = 1.2093375037721e+00; + P6 = 1.0236491595704e+00; + P7 = 9.9685258909811e-01; + P8 = 1.0004766563923e+00; + P9 = 9.9993617879146e-01; + P10 = 1.0000063122914e+00; + P11 = 9.9999966373260e-01; + + for i = 0:BP - 1 + P(i + 1) = eval(['P' num2str(i)]); + end + + Hv = ones(N, 1); + Hv(1:BP) = P; + Hv(end - BP + 1:end) = flip(P); + Hv = h * Hv; + H = spdiags(Hv, 0, N, N); + HI = spdiags(1 ./ Hv, 0, N, N); + %%%%%%%%%%%%%%%%%%%%%%%%% + + %%%% Q matrix %%%%%%%%%%% + + % interior stencil + 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); + + % Boundaries + Q0_0 = -5.0000000000000e-01; + Q0_1 = 6.7597560728423e-01; + Q0_2 = -2.6859785384416e-01; + Q0_3 = 1.4850302678903e-01; + Q0_4 = -8.7976689586154e-02; + Q0_5 = 4.1833336322613e-02; + Q0_6 = -2.2216684976993e-03; + Q0_7 = -1.5910034062022e-02; + Q0_8 = 1.1296706376589e-02; + Q0_9 = -3.1823678285130e-03; + Q0_10 = 2.4843594063649e-04; + Q0_11 = 3.1501105449828e-05; + Q0_12 = 0.0000000000000e+00; + Q0_13 = 0.0000000000000e+00; + Q0_14 = 0.0000000000000e+00; + Q0_15 = 0.0000000000000e+00; + Q0_16 = 0.0000000000000e+00; + Q0_17 = 0.0000000000000e+00; + Q1_0 = -6.7597560728423e-01; + Q1_1 = 0.0000000000000e+00; + Q1_2 = 9.5424013647146e-01; + Q1_3 = -4.3389334603464e-01; + Q1_4 = 2.4285669347653e-01; + Q1_5 = -1.1443465137214e-01; + Q1_6 = 8.5942765682435e-03; + Q1_7 = 4.0290424215772e-02; + Q1_8 = -2.9396383714543e-02; + Q1_9 = 8.5601827834256e-03; + Q1_10 = -7.8128092862319e-04; + Q1_11 = -6.0444181254875e-05; + Q1_12 = 0.0000000000000e+00; + Q1_13 = 0.0000000000000e+00; + Q1_14 = 0.0000000000000e+00; + Q1_15 = 0.0000000000000e+00; + Q1_16 = 0.0000000000000e+00; + Q1_17 = 0.0000000000000e+00; + Q2_0 = 2.6859785384416e-01; + Q2_1 = -9.5424013647146e-01; + Q2_2 = 0.0000000000000e+00; + Q2_3 = 9.7065114311923e-01; + Q2_4 = -4.3205328628292e-01; + Q2_5 = 1.9549970932735e-01; + Q2_6 = -2.4406885385172e-02; + Q2_7 = -5.5737279079895e-02; + Q2_8 = 4.3772338637753e-02; + Q2_9 = -1.3727655130726e-02; + Q2_10 = 1.6271304373071e-03; + Q2_11 = 1.7066984372933e-05; + Q2_12 = 0.0000000000000e+00; + Q2_13 = 0.0000000000000e+00; + Q2_14 = 0.0000000000000e+00; + Q2_15 = 0.0000000000000e+00; + Q2_16 = 0.0000000000000e+00; + Q2_17 = 0.0000000000000e+00; + Q3_0 = -1.4850302678903e-01; + Q3_1 = 4.3389334603464e-01; + Q3_2 = -9.7065114311923e-01; + Q3_3 = 0.0000000000000e+00; + Q3_4 = 9.5375878629204e-01; + Q3_5 = -3.6113954384951e-01; + Q3_6 = 6.9749289223875e-02; + Q3_7 = 6.5161366516465e-02; + Q3_8 = -6.0325702283960e-02; + Q3_9 = 2.1188913621662e-02; + Q3_10 = -3.2632650250470e-03; + Q3_11 = 1.3097937809499e-04; + Q3_12 = 0.0000000000000e+00; + Q3_13 = 0.0000000000000e+00; + Q3_14 = 0.0000000000000e+00; + Q3_15 = 0.0000000000000e+00; + Q3_16 = 0.0000000000000e+00; + Q3_17 = 0.0000000000000e+00; + Q4_0 = 8.7976689586154e-02; + Q4_1 = -2.4285669347653e-01; + Q4_2 = 4.3205328628292e-01; + Q4_3 = -9.5375878629204e-01; + Q4_4 = 0.0000000000000e+00; + Q4_5 = 8.8676146394834e-01; + Q4_6 = -2.1292503103800e-01; + Q4_7 = -4.6037018833218e-02; + Q4_8 = 7.4338719466734e-02; + Q4_9 = -3.1217656663809e-02; + Q4_10 = 6.1239492854797e-03; + Q4_11 = -4.5892226603067e-04; + Q4_12 = 0.0000000000000e+00; + Q4_13 = 0.0000000000000e+00; + Q4_14 = 0.0000000000000e+00; + Q4_15 = 0.0000000000000e+00; + Q4_16 = 0.0000000000000e+00; + Q4_17 = 0.0000000000000e+00; + Q5_0 = -4.1833336322613e-02; + Q5_1 = 1.1443465137214e-01; + Q5_2 = -1.9549970932735e-01; + Q5_3 = 3.6113954384951e-01; + Q5_4 = -8.8676146394834e-01; + Q5_5 = 0.0000000000000e+00; + Q5_6 = 7.7461223007026e-01; + Q5_7 = -1.0609547334165e-01; + Q5_8 = -4.4853791547749e-02; + Q5_9 = 3.2436468405486e-02; + Q5_10 = -8.4387621360184e-03; + Q5_11 = 8.5964292632428e-04; + Q5_12 = 0.0000000000000e+00; + Q5_13 = 0.0000000000000e+00; + Q5_14 = 0.0000000000000e+00; + Q5_15 = 0.0000000000000e+00; + Q5_16 = 0.0000000000000e+00; + Q5_17 = 0.0000000000000e+00; + Q6_0 = 2.2216684976993e-03; + Q6_1 = -8.5942765682435e-03; + Q6_2 = 2.4406885385172e-02; + Q6_3 = -6.9749289223875e-02; + Q6_4 = 2.1292503103800e-01; + Q6_5 = -7.7461223007026e-01; + Q6_6 = 0.0000000000000e+00; + Q6_7 = 7.4758103262966e-01; + Q6_8 = -1.5730779067906e-01; + Q6_9 = 2.6517620342970e-02; + Q6_10 = -4.3175367549700e-03; + Q6_11 = 1.1092605832824e-03; + Q6_12 = -1.8037518037522e-04; + Q6_13 = 0.0000000000000e+00; + Q6_14 = 0.0000000000000e+00; + Q6_15 = 0.0000000000000e+00; + Q6_16 = 0.0000000000000e+00; + Q6_17 = 0.0000000000000e+00; + Q7_0 = 1.5910034062022e-02; + Q7_1 = -4.0290424215772e-02; + Q7_2 = 5.5737279079895e-02; + Q7_3 = -6.5161366516465e-02; + Q7_4 = 4.6037018833218e-02; + Q7_5 = 1.0609547334165e-01; + Q7_6 = -7.4758103262966e-01; + Q7_7 = 0.0000000000000e+00; + Q7_8 = 8.0975719267918e-01; + Q7_9 = -2.3568822398349e-01; + Q7_10 = 6.9373143801571e-02; + Q7_11 = -1.6606121869177e-02; + Q7_12 = 2.5974025974031e-03; + Q7_13 = -1.8037518037522e-04; + Q7_14 = 0.0000000000000e+00; + Q7_15 = 0.0000000000000e+00; + Q7_16 = 0.0000000000000e+00; + Q7_17 = 0.0000000000000e+00; + Q8_0 = -1.1296706376589e-02; + Q8_1 = 2.9396383714543e-02; + Q8_2 = -4.3772338637753e-02; + Q8_3 = 6.0325702283960e-02; + Q8_4 = -7.4338719466734e-02; + Q8_5 = 4.4853791547749e-02; + Q8_6 = 1.5730779067906e-01; + Q8_7 = -8.0975719267918e-01; + Q8_8 = 0.0000000000000e+00; + Q8_9 = 8.4765775072084e-01; + Q8_10 = -2.6369594097148e-01; + Q8_11 = 7.8759594625702e-02; + Q8_12 = -1.7857142857146e-02; + Q8_13 = 2.5974025974031e-03; + Q8_14 = -1.8037518037522e-04; + Q8_15 = 0.0000000000000e+00; + Q8_16 = 0.0000000000000e+00; + Q8_17 = 0.0000000000000e+00; + Q9_0 = 3.1823678285130e-03; + Q9_1 = -8.5601827834256e-03; + Q9_2 = 1.3727655130726e-02; + Q9_3 = -2.1188913621662e-02; + Q9_4 = 3.1217656663809e-02; + Q9_5 = -3.2436468405486e-02; + Q9_6 = -2.6517620342970e-02; + Q9_7 = 2.3568822398349e-01; + Q9_8 = -8.4765775072084e-01; + Q9_9 = 0.0000000000000e+00; + Q9_10 = 8.5631774953989e-01; + Q9_11 = -2.6769768119702e-01; + Q9_12 = 7.9365079365093e-02; + Q9_13 = -1.7857142857146e-02; + Q9_14 = 2.5974025974031e-03; + Q9_15 = -1.8037518037522e-04; + Q9_16 = 0.0000000000000e+00; + Q9_17 = 0.0000000000000e+00; + Q10_0 = -2.4843594063649e-04; + Q10_1 = 7.8128092862319e-04; + Q10_2 = -1.6271304373071e-03; + Q10_3 = 3.2632650250470e-03; + Q10_4 = -6.1239492854797e-03; + Q10_5 = 8.4387621360184e-03; + Q10_6 = 4.3175367549700e-03; + Q10_7 = -6.9373143801571e-02; + Q10_8 = 2.6369594097148e-01; + Q10_9 = -8.5631774953989e-01; + Q10_10 = 0.0000000000000e+00; + Q10_11 = 8.5712580212095e-01; + Q10_12 = -2.6785714285718e-01; + Q10_13 = 7.9365079365093e-02; + Q10_14 = -1.7857142857146e-02; + Q10_15 = 2.5974025974031e-03; + Q10_16 = -1.8037518037522e-04; + Q10_17 = 0.0000000000000e+00; + Q11_0 = -3.1501105449828e-05; + Q11_1 = 6.0444181254875e-05; + Q11_2 = -1.7066984372933e-05; + Q11_3 = -1.3097937809499e-04; + Q11_4 = 4.5892226603067e-04; + Q11_5 = -8.5964292632428e-04; + Q11_6 = -1.1092605832824e-03; + Q11_7 = 1.6606121869177e-02; + Q11_8 = -7.8759594625702e-02; + Q11_9 = 2.6769768119702e-01; + Q11_10 = -8.5712580212095e-01; + Q11_11 = 0.0000000000000e+00; + Q11_12 = 8.5714285714289e-01; + Q11_13 = -2.6785714285718e-01; + Q11_14 = 7.9365079365093e-02; + Q11_15 = -1.7857142857146e-02; + Q11_16 = 2.5974025974031e-03; + Q11_17 = -1.8037518037522e-04; + + for i = 1:BP + + for j = 1:BP + Q(i, j) = eval(['Q' num2str(i - 1) '_' num2str(j - 1)]); + Q(N + 1 - i, N + 1 - j) = -eval(['Q' num2str(i - 1) '_' num2str(j - 1)]); + end + + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%% + + %%% Undivided difference operators %%%% + % Closed with zeros at the first boundary nodes. + m = N; + + DD_6 = (diag(ones(m - 3, 1), 3) - 6 * diag(ones(m - 2, 1), 2) + 15 * diag(ones(m - 1, 1), 1) - 20 * diag(ones(m, 1), 0) + 15 * diag(ones(m - 1, 1), -1) - 6 * diag(ones(m - 2, 1), -2) + diag(ones(m - 3, 1), -3)); + DD_6(1:9, 1:12) = [0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0; 0.70504538217624485042e1 -0.15481773512478815978e2 0.15439628908176388584e2 -0.11355022328644135767e2 0.62553652361133177256e1 -0.23565505827757565576e1 0.44789845784655348861e0 0 0 0 0 0; 0 0.84178325750049836493e0 -0.30776567832262484696e1 0.55480476621608890820e1 -0.66451562589578036671e1 0.54681670852508283306e1 -0.26873907470793164527e1 0.55220578435115281189e0 0 0 0 0; 0 0 0.36124410255434790612e0 -0.18841869963752696615e1 0.49068751071627186058e1 -0.79710602635488243814e1 0.75771246506939812399e1 -0.36661050678180486359e1 0.67610846733109492706e0 0 0 0; 0 0 0 0.31883568286286899671e0 -0.22216567686182446798e1 0.72341822309820977868e1 -0.12215760254444741754e2 0.10698736280837039597e2 -0.46222617037529123987e1 0.80792453213389245238e0 0 0; 0 0 0 0 0.45451605753051033631e0 -0.36739526096585069959e1 0.11306936594922967126e2 -0.16769946247690595362e2 0.13179014442979029716e2 -0.54150410306154005497e1 0.91847279253199572926e0 0; 0 0 0 0 0 0.77355004999207313207e0 -0.54143198515959566716e1 0.14230335022999000858e2 -0.19303948318184081752e2 0.14605034473743418451e2 -0.58729419174319386168e1 0.98229054047748459948e0; ]; + DD_6(m - 8:m, m - 11:m) = [0.98229054047748459948e0 -0.58729419174319386168e1 0.14605034473743418451e2 -0.19303948318184081752e2 0.14230335022999000858e2 -0.54143198515959566716e1 0.77355004999207313207e0 0 0 0 0 0; 0 0.91847279253199572926e0 -0.54150410306154005497e1 0.13179014442979029716e2 -0.16769946247690595362e2 0.11306936594922967126e2 -0.36739526096585069959e1 0.45451605753051033631e0 0 0 0 0; 0 0 0.80792453213389245238e0 -0.46222617037529123987e1 0.10698736280837039597e2 -0.12215760254444741754e2 0.72341822309820977868e1 -0.22216567686182446798e1 0.31883568286286899671e0 0 0 0; 0 0 0 0.67610846733109492706e0 -0.36661050678180486359e1 0.75771246506939812399e1 -0.79710602635488243814e1 0.49068751071627186058e1 -0.18841869963752696615e1 0.36124410255434790612e0 0 0; 0 0 0 0 0.55220578435115281189e0 -0.26873907470793164527e1 0.54681670852508283306e1 -0.66451562589578036671e1 0.55480476621608890820e1 -0.30776567832262484696e1 0.84178325750049836493e0 0; 0 0 0 0 0 0.44789845784655348861e0 -0.23565505827757565576e1 0.62553652361133177256e1 -0.11355022328644135767e2 0.15439628908176388584e2 -0.15481773512478815978e2 0.70504538217624485042e1; 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0; ]; + DD_6 = sparse(DD_6); + + DD_7 = (-diag(ones(m - 4, 1), -4) + 7 * diag(ones(m - 3, 1), -3) - 21 * diag(ones(m - 2, 1), -2) + 35 * diag(ones(m - 1, 1), -1) - 35 * diag(ones(m, 1), 0) + 21 * diag(ones(m - 1, 1), 1) - 7 * diag(ones(m - 2, 1), 2) + diag(ones(m - 3, 1), 3)); + DD_7(1:10, 1:13) = [0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0; -0.70504538217624585762e1 0.16323556769979337662e2 -0.18517285691402663507e2 0.16903069990805048996e2 -0.12900521495071139822e2 0.78247176680265960663e1 -0.31352892049258744203e1 0.55220578435115360075e0 0 0 0 0 0; 0 -0.77136636008199086135e0 0.31512297676528003033e1 -0.68105129732006589582e1 0.10585680229488408490e2 -0.12315008394369015993e2 0.94058676147776075845e1 -0.38654404904580696832e1 0.61955060619091911792e0 0 0 0 0; 0 0 -0.32268062071305431283e0 0.19678458985349116625e1 -0.63675477999083230026e1 0.13582054493165302513e2 -0.17679957518285956226e2 0.12831367737363170226e2 -0.47327592713176644894e1 0.72167708116161363042e0 0 0 0; 0 0 0 -0.28975994984220671445e0 0.24321233335964448997e1 -0.99133841479577475282e1 0.21377580445278298070e2 -0.24963717988619759059e2 0.16177915963135193396e2 -0.56554717249372471666e1 0.83471406934702410357e0 0 0; 0 0 0 0 -0.43028213606032104517e0 0.42103703770684731501e1 -0.15829711232892153976e2 0.29347405933458541883e2 -0.30751033700284402671e2 0.18952643607153901924e2 -0.64293095477239701048e1 0.92991669927993083983e0 0; 0 0 0 0 0 -0.76177813656089336989e0 0.63167064935286161169e1 -0.19922469032198601202e2 0.33781909556822143066e2 -0.34078413772067976385e2 0.20555296711011785159e2 -0.68760337833423921963e1 0.98478196280731881078e0; ]; + DD_7(m - 8:m, m - 12:m) = [-0.98478196280731881078e0 0.68760337833423921963e1 -0.20555296711011785159e2 0.34078413772067976385e2 -0.33781909556822143066e2 0.19922469032198601202e2 -0.63167064935286161169e1 0.76177813656089336989e0 0 0 0 0 0; 0 -0.92991669927993083983e0 0.64293095477239701048e1 -0.18952643607153901924e2 0.30751033700284402671e2 -0.29347405933458541883e2 0.15829711232892153976e2 -0.42103703770684731501e1 0.43028213606032104517e0 0 0 0 0; 0 0 -0.83471406934702410357e0 0.56554717249372471666e1 -0.16177915963135193396e2 0.24963717988619759059e2 -0.21377580445278298070e2 0.99133841479577475282e1 -0.24321233335964448997e1 0.28975994984220671445e0 0 0 0; 0 0 0 -0.72167708116161363042e0 0.47327592713176644894e1 -0.12831367737363170226e2 0.17679957518285956226e2 -0.13582054493165302513e2 0.63675477999083230026e1 -0.19678458985349116625e1 0.32268062071305431283e0 0 0; 0 0 0 0 -0.61955060619091911792e0 0.38654404904580696832e1 -0.94058676147776075845e1 0.12315008394369015993e2 -0.10585680229488408490e2 0.68105129732006589582e1 -0.31512297676528003033e1 0.77136636008199086135e0 0; 0 0 0 0 0 -0.55220578435115360075e0 0.31352892049258744203e1 -0.78247176680265960663e1 0.12900521495071139822e2 -0.16903069990805048996e2 0.18517285691402663507e2 -0.16323556769979337662e2 0.70504538217624585762e1; 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0; ]; + DD_7 = sparse(DD_7); + + DD_8 = (diag(ones(m - 4, 1), 4) - 8 * diag(ones(m - 3, 1), 3) + 28 * diag(ones(m - 2, 1), 2) - 56 * diag(ones(m - 1, 1), 1) + 70 * diag(ones(m, 1), 0) - 56 * diag(ones(m - 1, 1), -1) + 28 * diag(ones(m - 2, 1), -2) - 8 * diag(ones(m - 3, 1), -3) + diag(ones(m - 4, 1), -4)); + DD_8(1:10, 1:14) = [0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0.70504538217624673892e1 -0.17094923130061349891e2 0.21668515459055490895e2 -0.23713582964005737596e2 0.23486201724559577669e2 -0.20139726062395637234e2 0.12541156819703497681e2 -0.44176462748092288060e1 0.61955060619091989235e0 0 0 0 0 0; 0 0.71430915910501867497e0 -0.32169486987428931518e1 0.81290324137612243914e1 -0.15699214646211457760e2 0.23981482952559874949e2 -0.25082313639406953559e2 0.15461761961832278733e2 -0.49564048495273529434e1 0.66829534663026066516e0 0 0 0 0; 0 0 0.29213206789956748018e0 -0.20438756549159061384e1 0.79665959467484077329e1 -0.21271097908734749058e2 0.35359915036571912453e2 -0.34216980632968453935e2 0.18931037085270657958e2 -0.57734166492929090433e1 0.75569070942147255138e0 0 0 0; 0 0 0 0.26637215914130374043e0 -0.26313682263734604148e1 0.12983764630211125301e2 -0.34204128712445276912e2 0.49927435977239518119e2 -0.43141109235027182388e2 0.22621886899748988667e2 -0.66777125547761928286e1 0.85485906228117671607e0 0 0; 0 0 0 0 0.41007336094698233166e0 -0.47386249189431794318e1 0.21106281643856205301e2 -0.46955849493533667013e2 0.61502067400568805342e2 -0.50540382952410405130e2 0.25717238190895880419e2 -0.74393335942394467187e1 0.93853036285882489957e0 0; 0 0 0 0 0 0.75161513192516571838e0 -0.72190931354612755621e1 0.26563292042931468269e2 -0.54051055290915428906e2 0.68156827544135952769e2 -0.54814124562698093757e2 0.27504135133369568785e2 -0.78782557024585504862e1 0.98665883917119316896e0; ]; + DD_8(m - 9:m, m - 13:m) = [0.98665883917119316896e0 -0.78782557024585504862e1 0.27504135133369568785e2 -0.54814124562698093757e2 0.68156827544135952769e2 -0.54051055290915428906e2 0.26563292042931468269e2 -0.72190931354612755621e1 0.75161513192516571838e0 0 0 0 0 0; 0 0.93853036285882489957e0 -0.74393335942394467187e1 0.25717238190895880419e2 -0.50540382952410405130e2 0.61502067400568805342e2 -0.46955849493533667013e2 0.21106281643856205301e2 -0.47386249189431794318e1 0.41007336094698233166e0 0 0 0 0; 0 0 0.85485906228117671607e0 -0.66777125547761928286e1 0.22621886899748988667e2 -0.43141109235027182388e2 0.49927435977239518119e2 -0.34204128712445276912e2 0.12983764630211125301e2 -0.26313682263734604148e1 0.26637215914130374043e0 0 0 0; 0 0 0 0.75569070942147255138e0 -0.57734166492929090433e1 0.18931037085270657958e2 -0.34216980632968453935e2 0.35359915036571912453e2 -0.21271097908734749058e2 0.79665959467484077329e1 -0.20438756549159061384e1 0.29213206789956748018e0 0 0; 0 0 0 0 0.66829534663026066516e0 -0.49564048495273529434e1 0.15461761961832278733e2 -0.25082313639406953559e2 0.23981482952559874949e2 -0.15699214646211457760e2 0.81290324137612243914e1 -0.32169486987428931518e1 0.71430915910501867497e0 0; 0 0 0 0 0 0.61955060619091989235e0 -0.44176462748092288060e1 0.12541156819703497681e2 -0.20139726062395637234e2 0.23486201724559577669e2 -0.23713582964005737596e2 0.21668515459055490895e2 -0.17094923130061349891e2 0.70504538217624673892e1; 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0; ]; + DD_8 = sparse(DD_8); + + DD_9 = (-diag(ones(m - 5, 1), -5) + 9 * diag(ones(m - 4, 1), -4) - 36 * diag(ones(m - 3, 1), -3) + 84 * diag(ones(m - 2, 1), -2) - 126 * diag(ones(m - 1, 1), -1) + 126 * diag(ones(m, 1), 0) - 84 * diag(ones(m - 1, 1), 1) + 36 * diag(ones(m - 2, 1), 2) - 9 * diag(ones(m - 3, 1), 3) + diag(ones(m - 4, 1), 4)); + DD_9(1:11, 1:15) = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; -0.70504538217624752230e1 0.17809232289166388354e2 -0.24885464157798411697e2 0.31842615377766997368e2 -0.39185416370771078968e2 0.44121209014955561206e2 -0.37623470459110493043e2 0.19879408236641529627e2 -0.55759554557182790312e1 0.66829534663026140771e0 0 0 0 0 0; 0 -0.66695396914459764542e0 0.32764459415493351213e1 -0.94984942131335544919e1 0.22096883550779531311e2 -0.42252556942280502185e2 0.56435205688665645507e2 -0.46385285885496836199e2 0.22303821822873088245e2 -0.60146581196723459865e1 0.70559212586023632254e0 0 0 0 0; 0 0 -0.26728718140337933088e0 0.21137687176984662199e1 -0.96966416347745772182e1 0.31341597391980823551e2 -0.63647847065829442415e2 0.76988206424179021354e2 -0.56793111255811973873e2 0.25980374921818090695e2 -0.68012163847932529624e1 0.78215606693622397877e0 0 0 0; 0 0 0 -0.24708804973573377055e0 0.28212553166921198174e1 -0.16439370707786854975e2 0.51306193068667915368e2 -0.89869384759031132614e2 0.97067495778811160373e2 -0.67865660699246966000e2 0.30049706496492867728e2 -0.76937315605305904447e1 0.87058511566721451662e0 0 0; 0 0 0 0 -0.39286387086790423439e0 0.52598319320161875843e1 -0.27136647827815121102e2 0.70433774240300500520e2 -0.11070372132102384962e3 0.11371586164292341154e3 -0.77151714572687641258e2 0.33477001174077510234e2 -0.84467732657294240961e1 0.94525186880633042479e0 0; 0 0 0 0 0 -0.74268863896636254047e0 0.81214797773939350074e1 -0.34152804055197602060e2 0.81076582936373143359e2 -0.12268228957944471498e3 0.12333178026607071095e3 -0.82512405400108706356e2 0.35452150661063477188e2 -0.88799295525407385207e1 0.98812358535685795524e0; ]; + DD_9(m - 9:m, m - 14:m) = [-0.98812358535685795524e0 0.88799295525407385207e1 -0.35452150661063477188e2 0.82512405400108706356e2 -0.12333178026607071095e3 0.12268228957944471498e3 -0.81076582936373143359e2 0.34152804055197602060e2 -0.81214797773939350074e1 0.74268863896636254047e0 0 0 0 0 0; 0 -0.94525186880633042479e0 0.84467732657294240961e1 -0.33477001174077510234e2 0.77151714572687641258e2 -0.11371586164292341154e3 0.11070372132102384962e3 -0.70433774240300500520e2 0.27136647827815121102e2 -0.52598319320161875843e1 0.39286387086790423439e0 0 0 0 0; 0 0 -0.87058511566721451662e0 0.76937315605305904447e1 -0.30049706496492867728e2 0.67865660699246966000e2 -0.97067495778811160373e2 0.89869384759031132614e2 -0.51306193068667915368e2 0.16439370707786854975e2 -0.28212553166921198174e1 0.24708804973573377055e0 0 0 0; 0 0 0 -0.78215606693622397877e0 0.68012163847932529624e1 -0.25980374921818090695e2 0.56793111255811973873e2 -0.76988206424179021354e2 0.63647847065829442415e2 -0.31341597391980823551e2 0.96966416347745772182e1 -0.21137687176984662199e1 0.26728718140337933088e0 0 0; 0 0 0 0 -0.70559212586023632254e0 0.60146581196723459865e1 -0.22303821822873088245e2 0.46385285885496836199e2 -0.56435205688665645507e2 0.42252556942280502185e2 -0.22096883550779531311e2 0.94984942131335544919e1 -0.32764459415493351213e1 0.66695396914459764542e0 0; 0 0 0 0 0 -0.66829534663026140771e0 0.55759554557182790312e1 -0.19879408236641529627e2 0.37623470459110493043e2 -0.44121209014955561206e2 0.39185416370771078968e2 -0.31842615377766997368e2 0.24885464157798411697e2 -0.17809232289166388354e2 0.70504538217624752230e1; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; ]; + DD_9 = sparse(DD_9); + + DD_10 = (diag(ones(m - 5, 1), 5) - 10 * diag(ones(m - 4, 1), 4) + 45 * diag(ones(m - 3, 1), 3) - 120 * diag(ones(m - 2, 1), 2) + 210 * diag(ones(m - 1, 1), 1) - 252 * diag(ones(m, 1), 0) + 210 * diag(ones(m - 1, 1), -1) - 120 * diag(ones(m - 2, 1), -2) + 45 * diag(ones(m - 3, 1), -3) - 10 * diag(ones(m - 4, 1), -4) + diag(ones(m - 5, 1), -5)); + DD_10(1:11, 1:16) = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0.70504538217624822734e1 -0.18476186258311004476e2 0.28161910099347774980e2 -0.41341109590900593201e2 0.61282299921550671561e2 -0.86373765957236149764e2 0.94058676147776232608e2 -0.66264694122138432090e2 0.27879777278591395156e2 -0.66829534663026140771e1 0.70559212586023702812e0 0 0 0 0 0; 0 0.62689419647750179604e0 -0.33308831364980034247e1 0.10914786591113557343e2 -0.29883886064802801254e2 0.69173811658980680918e2 -0.11287041137733129101e3 0.11596321471374209050e3 -0.74346072742910294151e2 0.30073290598361729932e2 -0.70559212586023632254e1 0.73517682146919258207e0 0 0 0 0; 0 0 0.24665297575614270036e0 -0.21786018467637256442e1 0.11551532389532584562e2 -0.44092342567416599454e2 0.10607974510971573736e3 -0.15397641284835804271e3 0.14198277813952993468e3 -0.86601249739393635650e2 0.34006081923966264812e2 -0.78215606693622397877e1 0.80337713279357912969e0 0 0 0; 0 0 0 0.23087142251463535911e0 -0.30031734426541638399e1 0.20275063024062368195e2 -0.73294561526668450525e2 0.14978230793171855436e3 -0.19413499155762232075e3 0.16966415174811741500e3 -0.10016568832164289243e3 0.38468657802652952223e2 -0.87058511566721451662e1 0.88321407619404757308e0 0 0; 0 0 0 0 0.37796280007363075810e0 -0.57748488744870271355e1 0.33920809784768901377e2 -0.10061967748614357217e3 0.18450620220170641603e3 -0.22743172328584682309e3 0.19287928643171910314e3 -0.11159000391359170078e3 0.42233866328647120480e2 -0.94525186880633042479e1 0.95064470121725563376e0 0; 0 0 0 0 0 0.73474076934244039806e0 -0.90238664193265944527e1 0.42691005068997002575e2 -0.11582368990910449051e3 0.20447048263240785831e3 -0.24666356053214142190e3 0.20628101350027176589e3 -0.11817383553687825729e3 0.44399647762703692603e2 -0.98812358535685795524e1 0.98929851729658394154e0; ]; + DD_10(m - 10:m, m - 15:m) = [0.98929851729658394154e0 -0.98812358535685795524e1 0.44399647762703692603e2 -0.11817383553687825729e3 0.20628101350027176589e3 -0.24666356053214142190e3 0.20447048263240785831e3 -0.11582368990910449051e3 0.42691005068997002575e2 -0.90238664193265944527e1 0.73474076934244039806e0 0 0 0 0 0; 0 0.95064470121725563376e0 -0.94525186880633042479e1 0.42233866328647120480e2 -0.11159000391359170078e3 0.19287928643171910314e3 -0.22743172328584682309e3 0.18450620220170641603e3 -0.10061967748614357217e3 0.33920809784768901377e2 -0.57748488744870271355e1 0.37796280007363075810e0 0 0 0 0; 0 0 0.88321407619404757308e0 -0.87058511566721451662e1 0.38468657802652952223e2 -0.10016568832164289243e3 0.16966415174811741500e3 -0.19413499155762232075e3 0.14978230793171855436e3 -0.73294561526668450525e2 0.20275063024062368195e2 -0.30031734426541638399e1 0.23087142251463535911e0 0 0 0; 0 0 0 0.80337713279357912969e0 -0.78215606693622397877e1 0.34006081923966264812e2 -0.86601249739393635650e2 0.14198277813952993468e3 -0.15397641284835804271e3 0.10607974510971573736e3 -0.44092342567416599454e2 0.11551532389532584562e2 -0.21786018467637256442e1 0.24665297575614270036e0 0 0; 0 0 0 0 0.73517682146919258207e0 -0.70559212586023632254e1 0.30073290598361729932e2 -0.74346072742910294151e2 0.11596321471374209050e3 -0.11287041137733129101e3 0.69173811658980680918e2 -0.29883886064802801254e2 0.10914786591113557343e2 -0.33308831364980034247e1 0.62689419647750179604e0 0; 0 0 0 0 0 0.70559212586023702812e0 -0.66829534663026140771e1 0.27879777278591395156e2 -0.66264694122138432090e2 0.94058676147776232608e2 -0.86373765957236149764e2 0.61282299921550671561e2 -0.41341109590900593201e2 0.28161910099347774980e2 -0.18476186258311004476e2 0.70504538217624822734e1; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; ]; + DD_10 = sparse(DD_10); + + %[-1, 11, -55, 165, -330, 462, -462, 330, -165, 55, -11, 1, 0] + DD_11 = (-diag(ones(m - 6, 1), -6) + 11 * diag(ones(m - 5, 1), -5) - 55 * diag(ones(m - 4, 1), -4) + 165 * diag(ones(m - 3, 1), -3) - 330 * diag(ones(m - 2, 1), -2) + 462 * diag(ones(m - 1, 1), -1) - 462 * diag(ones(m, 1), 0) + 330 * diag(ones(m - 1, 1), 1) - 165 * diag(ones(m - 2, 1), 2) + 55 * diag(ones(m - 3, 1), 3) - 11 * diag(ones(m - 4, 1), 4) + diag(ones(m - 5, 1), 5)); + DD_11(1:12, 1:17) = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; -0.70504538217624886828e1 0.19103080454788523638e2 -0.31492793235845807035e2 0.52255896182014198048e2 -0.91166185986353555693e2 0.15554757761621697209e3 -0.20692908752510771174e3 0.18222790883588068825e3 -0.10222585002150178224e3 0.36756244064664377424e2 -0.77615133844626073094e1 0.73517682146919325041e0 0 0 0 0 0; 0 -0.59247568548574727441e0 0.33811178542850964393e1 -0.12374519230919224800e2 0.39160480492663455349e2 -0.10704747746062038309e3 0.20692908752510736686e3 -0.25511907237023259909e3 0.20445170004300330891e3 -0.11026873219399300975e3 0.38807566922312997740e2 -0.80869450361611184028e1 0.75926914003985711330e0 0 0 0 0; 0 0 -0.22922038452398770506e0 0.22391799149842653766e1 -0.13526028855965613027e2 0.59818136859120542509e2 -0.16669674231526758728e3 0.28229009022198974497e3 -0.31236211190696585630e3 0.23815343678333249804e3 -0.12468896705454297098e3 0.43018583681492318832e2 -0.88371484607293704266e1 0.82079151707601599340e0 0 0 0; 0 0 0 -0.21701391114437616116e0 0.31781915325611402199e1 -0.24486327517266714070e2 0.10078002209916911947e3 -0.23537219817841487113e3 0.35591415118897425470e3 -0.37326113384585831300e3 0.27545564288451795418e3 -0.14105174527639415815e3 0.47882181361696798414e2 -0.97153548381345233039e1 0.89358450029368883150e0 0 0; 0 0 0 0 -0.36488508570871329747e0 0.62843543720560487741e1 -0.41458767514717546128e2 0.13835205654344741174e3 -0.28993831774553865375e3 0.41695815935738584232e3 -0.42433443014978202692e3 0.30687251076237717714e3 -0.15485750987170610843e3 0.51988852784348173364e2 -0.10457091713389811971e2 0.95506826122820715686e0 0; 0 0 0 0 0 -0.72758579432539352184e0 0.99262530612592538979e1 -0.52177895084329669814e2 0.15925757362501867446e3 -0.32131075842235520591e3 0.45221652764225927349e3 -0.45381822970059788496e3 0.32497804772641520756e3 -0.16279870846324687288e3 0.54346797194627187538e2 -0.10882283690262423357e2 0.99026190553785350209e0; ]; + DD_11(m - 10:m, m - 16:m) = [-0.99026190553785350209e0 0.10882283690262423357e2 -0.54346797194627187538e2 0.16279870846324687288e3 -0.32497804772641520756e3 0.45381822970059788496e3 -0.45221652764225927349e3 0.32131075842235520591e3 -0.15925757362501867446e3 0.52177895084329669814e2 -0.99262530612592538979e1 0.72758579432539352184e0 0 0 0 0 0; 0 -0.95506826122820715686e0 0.10457091713389811971e2 -0.51988852784348173364e2 0.15485750987170610843e3 -0.30687251076237717714e3 0.42433443014978202692e3 -0.41695815935738584232e3 0.28993831774553865375e3 -0.13835205654344741174e3 0.41458767514717546128e2 -0.62843543720560487741e1 0.36488508570871329747e0 0 0 0 0; 0 0 -0.89358450029368883150e0 0.97153548381345233039e1 -0.47882181361696798414e2 0.14105174527639415815e3 -0.27545564288451795418e3 0.37326113384585831300e3 -0.35591415118897425470e3 0.23537219817841487113e3 -0.10078002209916911947e3 0.24486327517266714070e2 -0.31781915325611402199e1 0.21701391114437616116e0 0 0 0; 0 0 0 -0.82079151707601599340e0 0.88371484607293704266e1 -0.43018583681492318832e2 0.12468896705454297098e3 -0.23815343678333249804e3 0.31236211190696585630e3 -0.28229009022198974497e3 0.16669674231526758728e3 -0.59818136859120542509e2 0.13526028855965613027e2 -0.22391799149842653766e1 0.22922038452398770506e0 0 0; 0 0 0 0 -0.75926914003985711330e0 0.80869450361611184028e1 -0.38807566922312997740e2 0.11026873219399300975e3 -0.20445170004300330891e3 0.25511907237023259909e3 -0.20692908752510736686e3 0.10704747746062038309e3 -0.39160480492663455349e2 0.12374519230919224800e2 -0.33811178542850964393e1 0.59247568548574727441e0 0; 0 0 0 0 0 -0.73517682146919325041e0 0.77615133844626073094e1 -0.36756244064664377424e2 0.10222585002150178224e3 -0.18222790883588068825e3 0.20692908752510771174e3 -0.15554757761621697209e3 0.91166185986353555693e2 -0.52255896182014198048e2 0.31492793235845807035e2 -0.19103080454788523638e2 0.70504538217624886828e1; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; ]; + DD_11 = sparse(DD_11); + + %[1, -12, 66, -220, 495, -792, 924, -792, 495, -220, 66,-12, 1] + DD_12 = (diag(ones(m - 6, 1), 6) - 12 * diag(ones(m - 5, 1), 5) + 66 * diag(ones(m - 4, 1), 4) - 220 * diag(ones(m - 3, 1), 3) + 495 * diag(ones(m - 2, 1), 2) - 792 * diag(ones(m - 1, 1), 1) + 924 * diag(ones(m, 1), 0) - 792 * diag(ones(m - 1, 1), -1) + 495 * diag(ones(m - 2, 1), -2) - 220 * diag(ones(m - 3, 1), -3) + 66 * diag(ones(m - 4, 1), -4) - 12 * diag(ones(m - 5, 1), -5) + diag(ones(m - 6, 1), -6)); + DD_12(1:12, 1:18) = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0.70504538217624945581e1 -0.19695556140274287325e2 0.34873911090130932535e2 -0.64630415412933476707e2 0.13032666647901711965e3 -0.26259505507683757401e3 0.41385817505021542347e3 -0.43734698120611365179e3 0.30667755006450534672e3 -0.14702497625865750970e3 0.46569080306775643856e2 -0.88221218576303190049e1 0.75926914003985774602e0 0 0 0 0 0; 0 0.56252054413792412864e0 -0.34278021535209417538e1 0.13874841106233553089e2 -0.50022717612825328606e2 0.15842900977125025845e3 -0.35473557861446977176e3 0.51023814474046519819e3 -0.49068408010320794140e3 0.33080619658197902926e3 -0.15523026768925199096e3 0.48521670216966710417e2 -0.91112296804782853596e1 0.77929289272158630803e0 0 0 0 0; 0 0 0.21428192906429919385e0 -0.22961219278627835050e1 0.15615594467315483385e2 -0.78810282483468544556e2 0.25004511347290138092e3 -0.48392586895198241994e3 0.62472422381393171260e3 -0.57156824827999799529e3 0.37406690116362891293e3 -0.17207433472596927533e3 0.53022890764376222560e2 -0.98494982049121919208e1 0.83534896297519895228e0 0 0 0; 0 0 0 0.20501361895561890151e0 -0.33471539032596363119e1 0.29069145008244604383e2 -0.13437336279889215930e3 0.35305829726762230670e3 -0.61013854489538443663e3 0.74652226769171662600e3 -0.66109354292284309003e3 0.42315523582918247446e3 -0.19152872544678719366e3 0.58292129028807139823e2 -0.10723014003524265978e2 0.90225552616201164306e0 0 0; 0 0 0 0 0.35327850260839005288e0 -0.67888982569607603271e1 0.49750521017661055353e2 -0.18446940872459654898e3 0.43490747661830798063e3 -0.71478541604123287255e3 0.84866886029956405384e3 -0.73649402582970522515e3 0.46457252961511832529e3 -0.20795541113739269345e3 0.62742550280338871828e2 -0.11460819134738485882e2 0.95876279102790935799e0 0; 0 0 0 0 0 0.72108566182178027310e0 -0.10828639703191913343e2 0.62613474101195603777e2 -0.21234343150002489927e3 0.48196613763353280887e3 -0.77522833310101589741e3 0.90763645940119576992e3 -0.77994731454339649814e3 0.48839612538974061864e3 -0.21738718877850875015e3 0.65293702141574540142e2 -0.11883142866454242025e2 0.99106616353107873319e0; ]; + DD_12(m - 11:m, m - 17:m) = [0.99106616353107873319e0 -0.11883142866454242025e2 0.65293702141574540142e2 -0.21738718877850875015e3 0.48839612538974061864e3 -0.77994731454339649814e3 0.90763645940119576992e3 -0.77522833310101589741e3 0.48196613763353280887e3 -0.21234343150002489927e3 0.62613474101195603777e2 -0.10828639703191913343e2 0.72108566182178027310e0 0 0 0 0 0; 0 0.95876279102790935799e0 -0.11460819134738485882e2 0.62742550280338871828e2 -0.20795541113739269345e3 0.46457252961511832529e3 -0.73649402582970522515e3 0.84866886029956405384e3 -0.71478541604123287255e3 0.43490747661830798063e3 -0.18446940872459654898e3 0.49750521017661055353e2 -0.67888982569607603271e1 0.35327850260839005288e0 0 0 0 0; 0 0 0.90225552616201164306e0 -0.10723014003524265978e2 0.58292129028807139823e2 -0.19152872544678719366e3 0.42315523582918247446e3 -0.66109354292284309003e3 0.74652226769171662600e3 -0.61013854489538443663e3 0.35305829726762230670e3 -0.13437336279889215930e3 0.29069145008244604383e2 -0.33471539032596363119e1 0.20501361895561890151e0 0 0 0; 0 0 0 0.83534896297519895228e0 -0.98494982049121919208e1 0.53022890764376222560e2 -0.17207433472596927533e3 0.37406690116362891293e3 -0.57156824827999799529e3 0.62472422381393171260e3 -0.48392586895198241994e3 0.25004511347290138092e3 -0.78810282483468544556e2 0.15615594467315483385e2 -0.22961219278627835050e1 0.21428192906429919385e0 0 0; 0 0 0 0 0.77929289272158630803e0 -0.91112296804782853596e1 0.48521670216966710417e2 -0.15523026768925199096e3 0.33080619658197902926e3 -0.49068408010320794140e3 0.51023814474046519819e3 -0.35473557861446977176e3 0.15842900977125025845e3 -0.50022717612825328606e2 0.13874841106233553089e2 -0.34278021535209417538e1 0.56252054413792412864e0 0; 0 0 0 0 0 0.75926914003985774602e0 -0.88221218576303190049e1 0.46569080306775643856e2 -0.14702497625865750970e3 0.30667755006450534672e3 -0.43734698120611365179e3 0.41385817505021542347e3 -0.26259505507683757401e3 0.13032666647901711965e3 -0.64630415412933476707e2 0.34873911090130932535e2 -0.19695556140274287325e2 0.70504538217624945581e1; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; ]; + DD_12 = sparse(DD_12); + %%%%%%%%%%%%%%%%%%%%%%%%%%% + + %%%% Difference operators %% + D1 = H \ Q; + + % Helper functions for constructing D2(c) + % TODO: Consider changing sparse(diag(...)) to spdiags(....) + + % Minimal 13 point stencil width + function D2 = D2_fun_minimal(c) + % Here we add variable diffusion + C1 = sparse(diag(c)); + C2 = 1/2 * diag(ones(m - 1, 1), -1) + 1/2 * diag(ones(m, 1), 0); C2(1, 2) = 1/2; + C3 = 5/18 * diag(ones(m - 1, 1), -1) + 5/18 * diag(ones(m - 1, 1), 1) + 4/9 * diag(ones(m, 1), 0); C3(1, 3) = 5/18; C3(m, m - 2) = 5/18; + C4 = 5/28 * diag(ones(m - 2, 1), -2) + 9/28 * diag(ones(m - 1, 1), -1) + 5/28 * diag(ones(m - 1, 1), 1) + 9/28 * diag(ones(m, 1), 0); C4(2, 4) = 5/28; C4(1, 3) = 5/28; C4(1, 4) = 9/28; C4(m, m - 2) = 5/28; C4(m, m - 3) = 5/28; + C5 = 1/7 * diag(ones(m - 2, 1), -2) + 1/7 * diag(ones(m - 2, 1), 2) + 8/35 * diag(ones(m - 1, 1), -1) + 8/35 * diag(ones(m - 1, 1), 1) + 9/35 * diag(ones(m, 1), 0); C5(2, 4) = 1/7; C5(1, 3) = 1/7; C5(1, 4) = 1/7; C5(1, 5) = 8/35; C5(2, 5) = 1/7; C5(m - 1, m - 4) = 1/7; C5(m, m - 4) = 8/35; C5(m, m - 3) = 1/7; + C6 = 1/6 * diag(ones(m - 3, 1), -3) + 1/6 * diag(ones(m - 2, 1), -2) + 1/6 * diag(ones(m - 2, 1), 2) + 1/6 * diag(ones(m - 1, 1), -1) + 1/6 * diag(ones(m - 1, 1), 1) + 1/6 * diag(ones(m, 1), 0); + C6(1, 4) = 1/6; C6(1, 5) = 1/6; C6(1, 6) = 1/6; C6(2, 5) = 1/6; C6(2, 6) = 1/6; C6(3, 6) = 1/6; C6(m - 1, m - 5) = 1/6; C6(m, m - 4) = 1/6; C6(m, m - 5) = 1/6; + + C2 = sparse(diag(C2 * c)); + C3 = sparse(diag(C3 * c)); + C4 = sparse(diag(C4 * c)); + C5 = sparse(diag(C5 * c)); + C6 = sparse(diag(C6 * c)); + + % Remainder term added to wide second derivative operator + R = (1/30735936 / h) * transpose(DD_12) * C1 * DD_12 + (1/6403320 / h) * transpose(DD_11) * C2 * DD_11 + (1/1293600 / h) * transpose(DD_10) * C3 * DD_10 + (1/249480 / h) * transpose(DD_9) * C4 * DD_9 + (1/44352 / h) * transpose(DD_8) * C5 * DD_8 + (1/6468 / h) * transpose(DD_7) * C6 * DD_7; + D2 = D1 * C1 * D1 - H \ R; + end + + % Few additional grid point in interior stencil cmp. to minimal + function D2 = D2_fun_nonminimal(c) + % Here we add variable diffusion + C1 = sparse(diag(c)); + C2 = 1/2 * diag(ones(m - 1, 1), -1) + 1/2 * diag(ones(m, 1), 0); C2(1, 2) = 1/2; + C3 = 5/18 * diag(ones(m - 1, 1), -1) + 5/18 * diag(ones(m - 1, 1), 1) + 4/9 * diag(ones(m, 1), 0); C3(1, 3) = 5/18; C3(m, m - 2) = 5/18; + C4 = 5/28 * diag(ones(m - 2, 1), -2) + 9/28 * diag(ones(m - 1, 1), -1) + 5/28 * diag(ones(m - 1, 1), 1) + 9/28 * diag(ones(m, 1), 0); C4(2, 4) = 5/28; C4(1, 3) = 5/28; C4(1, 4) = 9/28; C4(m, m - 2) = 5/28; C4(m, m - 3) = 5/28; + C5 = 1/7 * diag(ones(m - 2, 1), -2) + 1/7 * diag(ones(m - 2, 1), 2) + 8/35 * diag(ones(m - 1, 1), -1) + 8/35 * diag(ones(m - 1, 1), 1) + 9/35 * diag(ones(m, 1), 0); C5(2, 4) = 1/7; C5(1, 3) = 1/7; C5(1, 4) = 1/7; C5(1, 5) = 8/35; C5(2, 5) = 1/7; C5(m - 1, m - 4) = 1/7; C5(m, m - 4) = 8/35; C5(m, m - 3) = 1/7; + + C2 = sparse(diag(C2 * c)); + C3 = sparse(diag(C3 * c)); + C4 = sparse(diag(C4 * c)); + C5 = sparse(diag(C5 * c)); + + % Remainder term added to wide second derivative operator + R = (1/30735936 / h) * transpose(DD_12) * C1 * DD_12 + (1/6403320 / h) * transpose(DD_11) * C2 * DD_11 + (1/1293600 / h) * transpose(DD_10) * C3 * DD_10 + (1/249480 / h) * transpose(DD_9) * C4 * DD_9 + (1/44352 / h) * transpose(DD_8) * C5 * DD_8; + D2 = D1 * C1 * D1 - H \ R; + end + + % Wide stencil + function D2 = D2_fun_wide(c) + % Here we add variable diffusion + C1 = sparse(diag(c)); + D2 = D1 * C1 * D1; + end + + switch options.stencil_width + case 'minimal' + D2 = @D2_fun_minimal; + case 'nonminimal' + D2 = @D2_fun_nonminimal; + case 'wide' + D2 = @D2_fun_wide; + otherwise + error('No option %s for stencil width', options.stencil_width) + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%% + + %%%% Artificial dissipation operator %%% + switch options.AD + case 'upwind' + % This is the choice that yield 11th order Upwind + DI = H \ (transpose(DD_6) * DD_6) * (-1/5544); + case 'op' + % This choice will preserve the order of the underlying + % Non-dissipative D1 SBP operator + DI = H \ (transpose(DD_7) * DD_7) * (-1 / (5 * 5544)); + % Notice that you can use any negative number instead of (-1/(5*5544)) + otherwise + error("Artificial dissipation options '%s' not implemented.", option.AD) + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%% +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/d2_noneq_variable_4.m Mon Feb 14 11:14:46 2022 +0100 @@ -0,0 +1,159 @@ +function [H, HI, D1, D2, DI] = d2_noneq_variable_4(N, h, options) + % N: Number of grid points + % h: grid spacing + % options: struct containing options for constructing the operator + % current options are: + % options.stencil_type ('minimal','nonminimal','wide') + % options.AD ('upwind', 'op') + + % BP: Number of boundary points + % order: Accuracy of interior stencil + BP = 4; + order = 4; + + %%%% Norm matrix %%%%%%%% + P = zeros(BP, 1); + P0 = 2.1259737557798e-01; + P1 = 1.0260290400758e+00; + P2 = 1.0775123588954e+00; + P3 = 9.8607273802835e-01; + + for i = 0:BP - 1 + P(i + 1) = eval(['P' num2str(i)]); + end + + Hv = ones(N, 1); + Hv(1:BP) = P; + Hv(end - BP + 1:end) = flip(P); + Hv = h * Hv; + H = spdiags(Hv, 0, N, N); + HI = spdiags(1 ./ Hv, 0, N, N); + %%%%%%%%%%%%%%%%%%%%%%%%% + + %%%% Q matrix %%%%%%%%%%% + d = [1/12, -2/3, 0, 2/3, -1/12]; + d = repmat(d, N, 1); + Q = spdiags(d, -order / 2:order / 2, N, N); + + % Boundaries + Q0_0 = -5.0000000000000e-01; + Q0_1 = 6.5605279837843e-01; + Q0_2 = -1.9875859409017e-01; + Q0_3 = 4.2705795711740e-02; + Q0_4 = 0.0000000000000e+00; + Q0_5 = 0.0000000000000e+00; + Q1_0 = -6.5605279837843e-01; + Q1_1 = 0.0000000000000e+00; + Q1_2 = 8.1236966439895e-01; + Q1_3 = -1.5631686602052e-01; + Q1_4 = 0.0000000000000e+00; + Q1_5 = 0.0000000000000e+00; + Q2_0 = 1.9875859409017e-01; + Q2_1 = -8.1236966439895e-01; + Q2_2 = 0.0000000000000e+00; + Q2_3 = 6.9694440364211e-01; + Q2_4 = -8.3333333333333e-02; + Q2_5 = 0.0000000000000e+00; + Q3_0 = -4.2705795711740e-02; + Q3_1 = 1.5631686602052e-01; + Q3_2 = -6.9694440364211e-01; + Q3_3 = 0.0000000000000e+00; + Q3_4 = 6.6666666666667e-01; + Q3_5 = -8.3333333333333e-02; + + for i = 1:BP + + for j = 1:BP + Q(i, j) = eval(['Q' num2str(i - 1) '_' num2str(j - 1)]); + Q(N + 1 - i, N + 1 - j) = -eval(['Q' num2str(i - 1) '_' num2str(j - 1)]); + end + + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%% + + %%% Undivided difference operators %%%% + % Closed with zeros at the first boundary nodes. + m = N; + + DD_2 = (diag(ones(m - 1, 1), -1) - 2 * diag(ones(m, 1), 0) + diag(ones(m - 1, 1), 1)); + DD_2(1:3, 1:4) = [0 0 0 0; 0.16138369498429727170e1 -0.26095138364100825853e1 0.99567688656710986834e0 0; 0 0.84859980956172494512e0 -0.17944203477786665350e1 0.94582053821694158989e0; ]; + DD_2(m - 2:m, m - 3:m) = [0.94582053821694158989e0 -0.17944203477786665350e1 0.84859980956172494512e0 0; 0 0.99567688656710986834e0 -0.26095138364100825853e1 0.16138369498429727170e1; 0 0 0 0; ]; + DD_2 = sparse(DD_2); + + DD_3 = (-diag(ones(m - 2, 1), -2) + 3 * diag(ones(m - 1, 1), -1) - 3 * diag(ones(m, 1), 0) + diag(ones(m - 1, 1), 1)); + DD_3(1:4, 1:5) = [0 0 0 0 0; 0 0 0 0 0; -0.17277463987989539852e1 0.37021976718569105700e1 -0.29870306597013296050e1 0.10125793866433730203e1 0; 0 -0.81738495424057284493e0 0.26916305216679998025e1 -0.28374616146508247697e1 0.96321604722339781208e0; ]; + DD_3(m - 2:m, m - 4:m) = [-0.96321604722339781208e0 0.28374616146508247697e1 -0.26916305216679998025e1 0.81738495424057284493e0 0; 0 -0.10125793866433730203e1 0.29870306597013296050e1 -0.37021976718569105700e1 0.17277463987989539852e1; 0 0 0 0 0; ]; + DD_3 = sparse(DD_3); + + DD_4 = (diag(ones(m - 2, 1), 2) - 4 * diag(ones(m - 1, 1), 1) + 6 * diag(ones(m, 1), 0) - 4 * diag(ones(m - 1, 1), -1) + diag(ones(m - 2, 1), -2)); + DD_4(1:4, 1:6) = [0 0 0 0 0 0; 0 0 0 0 0 0; 0.18176226052481525189e1 -0.47546882767009058782e1 0.59740613194026592100e1 -0.40503175465734920811e1 0.10133218986235862303e1 0; 0 0.79462567299107735362e0 -0.35888406955573330700e1 0.56749232293016495393e1 -0.38528641888935912483e1 0.97215598215819742539e0; ]; + DD_4(m - 3:m, m - 5:m) = [0.97215598215819742539e0 -0.38528641888935912483e1 0.56749232293016495393e1 -0.35888406955573330700e1 0.79462567299107735362e0 0; 0 0.10133218986235862303e1 -0.40503175465734920811e1 0.59740613194026592100e1 -0.47546882767009058782e1 0.18176226052481525189e1; 0 0 0 0 0 0; 0 0 0 0 0 0; ]; + DD_4 = sparse(DD_4); + %%%%%%%%%%%%%%%%%%%%%%%%%% + + %%%% Difference operators %%% + D1 = H \ Q; + + % Helper functions for constructing D2(c) + % TODO: Consider changing sparse(diag(...)) to spdiags(....) + + % Minimal 5 point stencil width + function D2 = D2_fun_minimal(c) + % Here we add variable diffusion + C1 = sparse(diag(c)); + C2 = 1/2 * diag(ones(m - 1, 1), -1) + 1/2 * diag(ones(m, 1), 0); C2(1, 2) = 1/2; + + C2 = sparse(diag(C2 * c)); + + % Remainder term added to wide second drivative opereator, to obtain a 5 + % point narrow stencil. + R = (1/144 / h) * transpose(DD_4) * C1 * DD_4 + (1/18 / h) * transpose(DD_3) * C2 * DD_3; + D2 = D1 * C1 * D1 - H \ R; + end + + % Few additional grid point in interior stencil cmp. to minimal + function D2 = D2_fun_nonminimal(c) + % Here we add variable diffusion + C1 = sparse(diag(c)); + + % Remainder term added to wide second derivative operator + R = (1/144 / h) * transpose(DD_4) * C1 * DD_4; + D2 = D1 * C1 * D1 - H \ R; + end + + % Wide stencil + function D2 = D2_fun_wide(c) + % Here we add variable diffusion + C1 = sparse(diag(c)); + D2 = D1 * C1 * D1; + end + + switch options.stencil_width + case 'minimal' + D2 = @D2_fun_minimal; + case 'nonminimal' + D2 = @D2_fun_nonminimal; + case 'wide' + D2 = @D2_fun_wide; + otherwise + error('No option %s for stencil width', options.stencil_width) + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%% + + %%%% Artificial dissipation operator %%% + switch options.AD + case 'upwind' + % This is the choice that yield 3rd order Upwind + DI = H \ (transpose(DD_2) * DD_2) * (-1/12); + case 'op' + % This choice will preserve the order of the underlying + % Non-dissipative D1 SBP operator + DI = H \ (transpose(DD_3) * DD_3) * (-1 / (5 * 12)); + otherwise + error("Artificial dissipation options '%s' not implemented.", option.AD) + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%% +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/d2_noneq_variable_6.m Mon Feb 14 11:14:46 2022 +0100 @@ -0,0 +1,202 @@ +function [H, HI, D1, D2, DI] = d2_noneq_variable_6(N, h, options) + % N: Number of grid points + % h: grid spacing + % options: struct containing options for constructing the operator + % current options are: + % options.stencil_type ('minimal','nonminimal','wide') + % options.AD ('upwind', 'op') + + % BP: Number of boundary points + % order: Accuracy of interior stencil + BP = 6; + order = 6; + + %%%% Norm matrix %%%%%%%% + P = zeros(BP, 1); + P0 = 1.3030223027124e-01; + P1 = 6.8851501587715e-01; + P2 = 9.5166202564389e-01; + P3 = 9.9103890475697e-01; + P4 = 1.0028757074552e+00; + P5 = 9.9950151111941e-01; + + for i = 0:BP - 1 + P(i + 1) = eval(['P' num2str(i)]); + end + + Hv = ones(N, 1); + Hv(1:BP) = P; + Hv(end - BP + 1:end) = flip(P); + Hv = h * Hv; + H = spdiags(Hv, 0, N, N); + HI = spdiags(1 ./ Hv, 0, N, N); + %%%%%%%%%%%%%%%%%%%%%%%%% + + %%%% Q matrix %%%%%%%%%%% + + % interior stencil + 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); + + % Boundaries + Q0_0 = -5.0000000000000e-01; + Q0_1 = 6.6042071945824e-01; + Q0_2 = -2.2104152954203e-01; + Q0_3 = 7.6243679810093e-02; + Q0_4 = -1.7298206716724e-02; + Q0_5 = 1.6753369904210e-03; + Q0_6 = 0.0000000000000e+00; + Q0_7 = 0.0000000000000e+00; + Q0_8 = 0.0000000000000e+00; + Q1_0 = -6.6042071945824e-01; + Q1_1 = 0.0000000000000e+00; + Q1_2 = 8.7352798702787e-01; + Q1_3 = -2.6581719253084e-01; + Q1_4 = 5.7458484948314e-02; + Q1_5 = -4.7485599871040e-03; + Q1_6 = 0.0000000000000e+00; + Q1_7 = 0.0000000000000e+00; + Q1_8 = 0.0000000000000e+00; + Q2_0 = 2.2104152954203e-01; + Q2_1 = -8.7352798702787e-01; + Q2_2 = 0.0000000000000e+00; + Q2_3 = 8.1707122038457e-01; + Q2_4 = -1.8881125503769e-01; + Q2_5 = 2.4226492138960e-02; + Q2_6 = 0.0000000000000e+00; + Q2_7 = 0.0000000000000e+00; + Q2_8 = 0.0000000000000e+00; + Q3_0 = -7.6243679810093e-02; + Q3_1 = 2.6581719253084e-01; + Q3_2 = -8.1707122038457e-01; + Q3_3 = 0.0000000000000e+00; + Q3_4 = 7.6798636652679e-01; + Q3_5 = -1.5715532552963e-01; + Q3_6 = 1.6666666666667e-02; + Q3_7 = 0.0000000000000e+00; + Q3_8 = 0.0000000000000e+00; + Q4_0 = 1.7298206716724e-02; + Q4_1 = -5.7458484948314e-02; + Q4_2 = 1.8881125503769e-01; + Q4_3 = -7.6798636652679e-01; + Q4_4 = 0.0000000000000e+00; + Q4_5 = 7.5266872305402e-01; + Q4_6 = -1.5000000000000e-01; + Q4_7 = 1.6666666666667e-02; + Q4_8 = 0.0000000000000e+00; + Q5_0 = -1.6753369904210e-03; + Q5_1 = 4.7485599871040e-03; + Q5_2 = -2.4226492138960e-02; + Q5_3 = 1.5715532552963e-01; + Q5_4 = -7.5266872305402e-01; + Q5_5 = 0.0000000000000e+00; + Q5_6 = 7.5000000000000e-01; + Q5_7 = -1.5000000000000e-01; + Q5_8 = 1.6666666666667e-02; + + for i = 1:BP + + for j = 1:BP + Q(i, j) = eval(['Q' num2str(i - 1) '_' num2str(j - 1)]); + Q(N + 1 - i, N + 1 - j) = -eval(['Q' num2str(i - 1) '_' num2str(j - 1)]); + end + + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%% + + %%% Undivided difference operators %%%% + % Closed with zeros at the first boundary nodes. + m = N; + + DD_3 = (-diag(ones(m - 2, 1), -2) + 3 * diag(ones(m - 1, 1), -1) - 3 * diag(ones(m, 1), 0) + diag(ones(m - 1, 1), 1)); + DD_3(1:5, 1:6) = [0 0 0 0 0 0; 0 0 0 0 0 0; -0.46757024540266021836e1 0.88373748766984018738e1 -0.56477423503490435435e1 0.14860699276772438533e1 0 0; 0 -0.13802450758054908946e1 0.36701915175801340778e1 -0.33643068661005748879e1 0.10743604243259317047e1 0; 0 0 -0.10409288946349185618e1 0.30665535320781497878e1 -0.30329117010471766032e1 0.10072870636039453772e1; ]; + DD_3(m - 3:m, m - 5:m) = [-0.10072870636039453772e1 0.30329117010471766032e1 -0.30665535320781497878e1 0.10409288946349185618e1 0 0; 0 -0.10743604243259317047e1 0.33643068661005748879e1 -0.36701915175801340778e1 0.13802450758054908946e1 0; 0 0 -0.14860699276772438533e1 0.56477423503490435435e1 -0.88373748766984018738e1 0.46757024540266021836e1; 0 0 0 0 0 0; ]; + DD_3 = sparse(DD_3); + + DD_4 = (diag(ones(m - 2, 1), 2) - 4 * diag(ones(m - 1, 1), 1) + 6 * diag(ones(m, 1), 0) - 4 * diag(ones(m - 1, 1), -1) + diag(ones(m - 2, 1), -2)); + DD_4(1:5, 1:7) = [0 0 0 0 0 0 0; 0 0 0 0 0 0 0; 0.57302111593550648941e1 -0.12521994384708052700e2 0.11419402572582197931e2 -0.59442797107089754133e1 0.13166603634797652881e1 0 0; 0 0.14441513881249918393e1 -0.49292485821432017638e1 0.67286137322011497757e1 -0.42974416973037268190e1 0.10539251591207869677e1 0; 0 0 0.10466075357769140419e1 -0.40887380427708663837e1 0.60658234020943532065e1 -0.40291482544157815088e1 0.10054553593153806442e1; ]; + DD_4(m - 4:m, m - 6:m) = [0.10054553593153806442e1 -0.40291482544157815088e1 0.60658234020943532065e1 -0.40887380427708663837e1 0.10466075357769140419e1 0 0; 0 0.10539251591207869677e1 -0.42974416973037268190e1 0.67286137322011497757e1 -0.49292485821432017638e1 0.14441513881249918393e1 0; 0 0 0.13166603634797652881e1 -0.59442797107089754133e1 0.11419402572582197931e2 -0.12521994384708052700e2 0.57302111593550648941e1; 0 0 0 0 0 0 0; 0 0 0 0 0 0 0; ]; + DD_4 = sparse(DD_4); + + DD_5 = (-diag(ones(m - 3, 1), -3) + 5 * diag(ones(m - 2, 1), -2) - 10 * diag(ones(m - 1, 1), -1) + 10 * diag(ones(m, 1), 0) - 5 * diag(ones(m - 1, 1), 1) + diag(ones(m - 2, 1), 2)); + DD_5(1:6, 1:8) = [0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0; -0.67194556014531368457e1 0.16377214352871472626e2 -0.19171027475746103125e2 0.14860699276772438533e2 -0.65833018173988264407e1 0.12358712649541552519e1 0 0; 0 -0.14971527633959360324e1 0.61951742553920293904e1 -0.11214356220335249626e2 0.10743604243259317047e2 -0.52696257956039348385e1 0.10423562806837740594e1 0; 0 0 -0.10511702536596915242e1 0.51109225534635829797e1 -0.10109705670157255344e2 0.10072870636039453772e2 -0.50272767965769032212e1 0.10043595308908133377e1; ]; + DD_5(m - 4:m, m - 7:m) = [-0.10043595308908133377e1 0.50272767965769032212e1 -0.10072870636039453772e2 0.10109705670157255344e2 -0.51109225534635829797e1 0.10511702536596915242e1 0 0; 0 -0.10423562806837740594e1 0.52696257956039348385e1 -0.10743604243259317047e2 0.11214356220335249626e2 -0.61951742553920293904e1 0.14971527633959360324e1 0; 0 0 -0.12358712649541552519e1 0.65833018173988264407e1 -0.14860699276772438533e2 0.19171027475746103125e2 -0.16377214352871472626e2 0.67194556014531368457e1; 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0; ]; + DD_5 = sparse(DD_5); + + DD_6 = (diag(ones(m - 3, 1), 3) - 6 * diag(ones(m - 2, 1), 2) + 15 * diag(ones(m - 1, 1), 1) - 20 * diag(ones(m, 1), 0) + 15 * diag(ones(m - 1, 1), -1) - 6 * diag(ones(m - 2, 1), -2) + diag(ones(m - 3, 1), -3)); + DD_6(1:6, 1:9) = [0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0; 0.76591061528436941127e1 -0.20373923615000091397e2 0.28913418478606999359e2 -0.29721398553544877066e2 0.19749905452196479322e2 -0.74152275897249315116e1 0.11881196746227271813e1 0 0; 0 0.15426631885693469226e1 -0.74666187707188589528e1 0.16821534330502874439e2 -0.21487208486518634095e2 0.15808877386811804515e2 -0.62541376841026443562e1 0.10348900354561115264e1 0; 0 0 0.10549863219420430611e1 -0.61331070641562995756e1 0.15164558505235883016e2 -0.20145741272078907544e2 0.15081830389730709664e2 -0.60261571853448800265e1 0.10036303046714514054e1; ]; + DD_6(m - 5:m, m - 8:m) = [0.10036303046714514054e1 -0.60261571853448800265e1 0.15081830389730709664e2 -0.20145741272078907544e2 0.15164558505235883016e2 -0.61331070641562995756e1 0.10549863219420430611e1 0 0; 0 0.10348900354561115264e1 -0.62541376841026443562e1 0.15808877386811804515e2 -0.21487208486518634095e2 0.16821534330502874439e2 -0.74666187707188589528e1 0.15426631885693469226e1 0; 0 0 0.11881196746227271813e1 -0.74152275897249315116e1 0.19749905452196479322e2 -0.29721398553544877066e2 0.28913418478606999359e2 -0.20373923615000091397e2 0.76591061528436941127e1; 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0; ]; + DD_6 = sparse(DD_6); + + %%%% Difference operators %%% + D1 = H \ Q; + + % Helper functions for constructing D2(c) + % TODO: Consider changing sparse(diag(...)) to spdiags(....) + + % Minimal 7 point stencil width + function D2 = D2_fun_minimal(c) + % Here we add variable diffusion + C1 = sparse(diag(c)); + C2 = 1/2 * diag(ones(m - 1, 1), -1) + 1/2 * diag(ones(m, 1), 0); C2(1, 2) = 1/2; + C3 = 1/3 * diag(ones(m - 1, 1), -1) + 1/3 * diag(ones(m - 1, 1), 1) + 1/3 * diag(ones(m, 1), 0); C3(1, 3) = 1/3; C3(m, m - 2) = 1/3; + + C2 = sparse(diag(C2 * c)); + C3 = sparse(diag(C3 * c)); + + % Remainder term added to wide second derivative operator + R = (1/3600 / h) * transpose(DD_6) * C1 * DD_6 + (1/600 / h) * transpose(DD_5) * C2 * DD_5 + (1/80 / h) * transpose(DD_4) * C3 * DD_4; + D2 = D1 * C1 * D1 - H \ R; + end + + % Few additional grid point in interior stencil cmp. to minimal + function D2 = D2_fun_nonminimal(c) + % Here we add variable diffusion + C1 = sparse(diag(c)); + C2 = 1/2 * diag(ones(m - 1, 1), -1) + 1/2 * diag(ones(m, 1), 0); C2(1, 2) = 1/2; + + C2 = sparse(diag(C2 * c)); + + % Remainder term added to wide second derivative operator + R = (1/3600 / h) * transpose(DD_6) * C1 * DD_6 + (1/600 / h) * transpose(DD_5) * C2 * DD_5; + D2 = D1 * C1 * D1 - H \ R; + end + + % Wide stencil + function D2 = D2_fun_wide(c) + % Here we add variable diffusion + C1 = sparse(diag(c)); + D2 = D1 * C1 * D1; + end + + switch options.stencil_width + case 'minimal' + D2 = @D2_fun_minimal; + case 'nonminimal' + D2 = @D2_fun_nonminimal; + case 'wide' + D2 = @D2_fun_wide; + otherwise + error('No option %s for stencil width', options.stencil_width) + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%% + + %%%% Artificial dissipation operator %%% + switch options.AD + case 'upwind' + % This is the choice that yield 3rd order Upwind + DI = H \ (transpose(DD_3) * DD_3) * (-1/60); + case 'op' + % This choice will preserve the order of the underlying + % Non-dissipative D1 SBP operator + DI = H \ (transpose(DD_4) * DD_4) * (-1 / (5 * 60)); + % Notice that you can use any negative number instead of (-1/(5*60)) + otherwise + error("Artificial dissipation options '%s' not implemented.", option.AD) + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%% +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/d2_noneq_variable_8.m Mon Feb 14 11:14:46 2022 +0100 @@ -0,0 +1,258 @@ +function [H, HI, D1, D2, DI] = d2_noneq_variable_8(N, h, options) + % N: Number of grid points + % h: grid spacing + % options: struct containing options for constructing the operator + % current options are: + % options.stencil_type ('minimal','nonminimal','wide') + % options.AD ('upwind', 'op') + + % BP: Number of boundary points + % order: Accuracy of interior stencil + BP = 8; + order = 8; + + %%%% Norm matrix %%%%%%%% + P = zeros(BP, 1); + P0 = 1.0758368078310e-01; + P1 = 6.1909685107891e-01; + P2 = 9.6971176519117e-01; + P3 = 1.1023441350947e+00; + P4 = 1.0244688965833e+00; + P5 = 9.9533550116831e-01; + P6 = 1.0008236941028e+00; + P7 = 9.9992060631812e-01; + + for i = 0:BP - 1 + P(i + 1) = eval(['P' num2str(i)]); + end + + Hv = ones(N, 1); + Hv(1:BP) = P; + Hv(end - BP + 1:end) = flip(P); + Hv = h * Hv; + H = spdiags(Hv, 0, N, N); + HI = spdiags(1 ./ Hv, 0, N, N); + %%%%%%%%%%%%%%%%%%%%%%%%% + + %%%% Q matrix %%%%%%%%%%% + + % interior stencil + 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); + + % Boundaries + Q0_0 = -5.0000000000000e-01; + Q0_1 = 6.7284756079369e-01; + Q0_2 = -2.5969732837062e-01; + Q0_3 = 1.3519390385721e-01; + Q0_4 = -6.9678474730984e-02; + Q0_5 = 2.6434024071371e-02; + Q0_6 = -5.5992311465618e-03; + Q0_7 = 4.9954552590464e-04; + Q0_8 = 0.0000000000000e+00; + Q0_9 = 0.0000000000000e+00; + Q0_10 = 0.0000000000000e+00; + Q0_11 = 0.0000000000000e+00; + Q1_0 = -6.7284756079369e-01; + Q1_1 = 0.0000000000000e+00; + Q1_2 = 9.4074021172233e-01; + Q1_3 = -4.0511642426516e-01; + Q1_4 = 1.9369192209331e-01; + Q1_5 = -6.8638079843479e-02; + Q1_6 = 1.3146457241484e-02; + Q1_7 = -9.7652615479254e-04; + Q1_8 = 0.0000000000000e+00; + Q1_9 = 0.0000000000000e+00; + Q1_10 = 0.0000000000000e+00; + Q1_11 = 0.0000000000000e+00; + Q2_0 = 2.5969732837062e-01; + Q2_1 = -9.4074021172233e-01; + Q2_2 = 0.0000000000000e+00; + Q2_3 = 9.4316393361096e-01; + Q2_4 = -3.5728039257451e-01; + Q2_5 = 1.1266686855013e-01; + Q2_6 = -1.8334941452280e-02; + Q2_7 = 8.2741521740941e-04; + Q2_8 = 0.0000000000000e+00; + Q2_9 = 0.0000000000000e+00; + Q2_10 = 0.0000000000000e+00; + Q2_11 = 0.0000000000000e+00; + Q3_0 = -1.3519390385721e-01; + Q3_1 = 4.0511642426516e-01; + Q3_2 = -9.4316393361096e-01; + Q3_3 = 0.0000000000000e+00; + Q3_4 = 8.7694387866575e-01; + Q3_5 = -2.4698058719506e-01; + Q3_6 = 4.7291642094198e-02; + Q3_7 = -4.0135203618880e-03; + Q3_8 = 0.0000000000000e+00; + Q3_9 = 0.0000000000000e+00; + Q3_10 = 0.0000000000000e+00; + Q3_11 = 0.0000000000000e+00; + Q4_0 = 6.9678474730984e-02; + Q4_1 = -1.9369192209331e-01; + Q4_2 = 3.5728039257451e-01; + Q4_3 = -8.7694387866575e-01; + Q4_4 = 0.0000000000000e+00; + Q4_5 = 8.1123946853807e-01; + Q4_6 = -2.0267150541446e-01; + Q4_7 = 3.8680398901392e-02; + Q4_8 = -3.5714285714286e-03; + Q4_9 = 0.0000000000000e+00; + Q4_10 = 0.0000000000000e+00; + Q4_11 = 0.0000000000000e+00; + Q5_0 = -2.6434024071371e-02; + Q5_1 = 6.8638079843479e-02; + Q5_2 = -1.1266686855013e-01; + Q5_3 = 2.4698058719506e-01; + Q5_4 = -8.1123946853807e-01; + Q5_5 = 0.0000000000000e+00; + Q5_6 = 8.0108544742793e-01; + Q5_7 = -2.0088756283071e-01; + Q5_8 = 3.8095238095238e-02; + Q5_9 = -3.5714285714286e-03; + Q5_10 = 0.0000000000000e+00; + Q5_11 = 0.0000000000000e+00; + Q6_0 = 5.5992311465618e-03; + Q6_1 = -1.3146457241484e-02; + Q6_2 = 1.8334941452280e-02; + Q6_3 = -4.7291642094198e-02; + Q6_4 = 2.0267150541446e-01; + Q6_5 = -8.0108544742793e-01; + Q6_6 = 0.0000000000000e+00; + Q6_7 = 8.0039405922650e-01; + Q6_8 = -2.0000000000000e-01; + Q6_9 = 3.8095238095238e-02; + Q6_10 = -3.5714285714286e-03; + Q6_11 = 0.0000000000000e+00; + Q7_0 = -4.9954552590464e-04; + Q7_1 = 9.7652615479254e-04; + Q7_2 = -8.2741521740941e-04; + Q7_3 = 4.0135203618880e-03; + Q7_4 = -3.8680398901392e-02; + Q7_5 = 2.0088756283071e-01; + Q7_6 = -8.0039405922650e-01; + Q7_7 = 0.0000000000000e+00; + Q7_8 = 8.0000000000000e-01; + Q7_9 = -2.0000000000000e-01; + Q7_10 = 3.8095238095238e-02; + Q7_11 = -3.5714285714286e-03; + + for i = 1:BP + + for j = 1:BP + Q(i, j) = eval(['Q' num2str(i - 1) '_' num2str(j - 1)]); + Q(N + 1 - i, N + 1 - j) = -eval(['Q' num2str(i - 1) '_' num2str(j - 1)]); + end + + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%% + + %%% Undivided difference operators %%%% + % Closed with zeros at the first boundary nodes. + m = N; + + DD_4 = (diag(ones(m - 2, 1), 2) - 4 * diag(ones(m - 1, 1), 1) + 6 * diag(ones(m, 1), 0) - 4 * diag(ones(m - 1, 1), -1) + diag(ones(m - 2, 1), -2)); + DD_4(1:6, 1:8) = [0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0; 0.70921010190504348684e1 -0.14196080536361841322e2 0.11072881931325435634e2 -0.50473576941871051066e1 0.10784552801730759259e1 0 0 0; 0 0.13740993382151221352e1 -0.42105600869792757010e1 0.54761010136211975317e1 -0.35797005751940657417e1 0.94006031033702177578e0 0 0; 0 0 0.82467928104463767301e0 -0.33274694995849432461e1 0.52587584638857303123e1 -0.37020511582893568152e1 0.94608291294393207601e0 0; 0 0 0 0.86436129166612654748e0 -0.37325441295306179390e1 0.57924699560798105338e1 -0.39066885960487908497e1 0.98240147783347170744e0; ]; + DD_4(m - 5:m, m - 7:m) = [0.98240147783347170744e0 -0.39066885960487908497e1 0.57924699560798105338e1 -0.37325441295306179390e1 0.86436129166612654748e0 0 0 0; 0 0.94608291294393207601e0 -0.37020511582893568152e1 0.52587584638857303123e1 -0.33274694995849432461e1 0.82467928104463767301e0 0 0; 0 0 0.94006031033702177578e0 -0.35797005751940657417e1 0.54761010136211975317e1 -0.42105600869792757010e1 0.13740993382151221352e1 0; 0 0 0 0.10784552801730759259e1 -0.50473576941871051066e1 0.11072881931325435634e2 -0.14196080536361841322e2 0.70921010190504348684e1; 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0; ]; + DD_4 = sparse(DD_4); + + DD_5 = (-diag(ones(m - 3, 1), -3) + 5 * diag(ones(m - 2, 1), -2) - 10 * diag(ones(m - 1, 1), -1) + 10 * diag(ones(m, 1), 0) - 5 * diag(ones(m - 1, 1), 1) + diag(ones(m - 2, 1), 2)); + DD_5(1:7, 1:9) = [0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0; -0.82098088052411907132e1 0.18024024120655590699e2 -0.17692096674769448317e2 0.12181944917152947702e2 -0.53922764008653796295e1 0.10882128430674802600e1 0 0 0; 0 -0.13913240333052950751e1 0.50983574122643719988e1 -0.89139255753021486176e1 0.89492514379851643542e1 -0.47003015516851088789e1 0.95794231004301621873e0 0 0; 0 0 -0.80388595981635823711e0 0.40861386922975635215e1 -0.87645974398095505205e1 0.92551278957233920380e1 -0.47304145647196603800e1 0.95763137632461357808e0 0; 0 0 0 -0.85214912336218661144e0 0.46656801619132724238e1 -0.96541165934663508896e1 0.97667214901219771242e1 -0.49120073891673585372e1 0.98587145396064649030e0; ]; + DD_5(m - 5:m, m - 8:m) = [-0.98587145396064649030e0 0.49120073891673585372e1 -0.97667214901219771242e1 0.96541165934663508896e1 -0.46656801619132724238e1 0.85214912336218661144e0 0 0 0; 0 -0.95763137632461357808e0 0.47304145647196603800e1 -0.92551278957233920380e1 0.87645974398095505205e1 -0.40861386922975635215e1 0.80388595981635823711e0 0 0; 0 0 -0.95794231004301621873e0 0.47003015516851088789e1 -0.89492514379851643542e1 0.89139255753021486176e1 -0.50983574122643719988e1 0.13913240333052950751e1 0; 0 0 0 -0.10882128430674802600e1 0.53922764008653796295e1 -0.12181944917152947702e2 0.17692096674769448317e2 -0.18024024120655590699e2 0.82098088052411907132e1; 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0; ]; + DD_5 = sparse(DD_5); + + DD_6 = (diag(ones(m - 3, 1), 3) - 6 * diag(ones(m - 2, 1), 2) + 15 * diag(ones(m - 1, 1), 1) - 20 * diag(ones(m, 1), 0) + 15 * diag(ones(m - 1, 1), -1) - 6 * diag(ones(m - 2, 1), -2) + diag(ones(m - 3, 1), -3)); + DD_6(1:7, 1:10) = [0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0; 0.92604272237009487731e1 -0.21899951980342041757e2 0.25706973995952182034e2 -0.23795532642767976509e2 0.16176829202596138889e2 -0.65292770584048815597e1 0.10805312592656301300e1 0 0 0; 0 0.14058275749850312569e1 -0.59637699688344396810e1 0.13135580487711615439e2 -0.17898502875970328708e2 0.14100904655055326637e2 -0.57476538602580973124e1 0.96761398731089236944e0 0 0; 0 0 0.78692381135906040550e0 -0.48340889923923043812e1 0.13146896159714325781e2 -0.18510255791446784076e2 0.14191243694158981140e2 -0.57457882579476814685e1 0.96506937655440259938e0 0; 0 0 0 0.84209241882516286974e0 -0.55988161942959269085e1 0.14481174890199526334e2 -0.19533442980243954248e2 0.14736022167502075612e2 -0.59152287237638789418e1 0.98819842177699528293e0; ]; + DD_6(m - 6:m, m - 9:m) = [0.98819842177699528293e0 -0.59152287237638789418e1 0.14736022167502075612e2 -0.19533442980243954248e2 0.14481174890199526334e2 -0.55988161942959269085e1 0.84209241882516286974e0 0 0 0; 0 0.96506937655440259938e0 -0.57457882579476814685e1 0.14191243694158981140e2 -0.18510255791446784076e2 0.13146896159714325781e2 -0.48340889923923043812e1 0.78692381135906040550e0 0 0; 0 0 0.96761398731089236944e0 -0.57476538602580973124e1 0.14100904655055326637e2 -0.17898502875970328708e2 0.13135580487711615439e2 -0.59637699688344396810e1 0.14058275749850312569e1 0; 0 0 0 0.10805312592656301300e1 -0.65292770584048815597e1 0.16176829202596138889e2 -0.23795532642767976509e2 0.25706973995952182034e2 -0.21899951980342041757e2 0.92604272237009487731e1; 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0; ]; + DD_6 = sparse(DD_6); + + DD_7 = (-diag(ones(m - 4, 1), -4) + 7 * diag(ones(m - 3, 1), -3) - 21 * diag(ones(m - 2, 1), -2) + 35 * diag(ones(m - 1, 1), -1) - 35 * diag(ones(m, 1), 0) + 21 * diag(ones(m - 1, 1), 1) - 7 * diag(ones(m - 2, 1), 2) + diag(ones(m - 3, 1), 3)); + DD_7(1:8, 1:11) = [0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0; -0.10257962606384161889e2 0.25816283570514673553e2 -0.35082323899232605336e2 0.40909341259657284245e2 -0.37745934806057657407e2 0.22852469704417085459e2 -0.75637188148594109098e1 0.10718455919447922848e1 0 0 0; 0 -0.14183700945143243060e1 0.68109221539151178574e1 -0.18129991367652287552e2 0.31322380032948075240e2 -0.32902110861795762152e2 0.20116788510903340593e2 -0.67732979111762465861e1 0.97367953737208690559e0 0 0; 0 0 -0.77264857229409862775e0 0.55732122985135573041e1 -0.18405654623600056093e2 0.32392947635031872133e2 -0.33112901953037622660e2 0.20110258902816885140e2 -0.67554856358808181957e1 0.97027194845028099974e0 0; 0 0 0 -0.83355973075426177031e0 0.65319522266785813933e1 -0.20273644846279336868e2 0.34183525215426919935e2 -0.34384051724171509761e2 0.20703300533173576296e2 -0.69173889524389669805e1 0.98986727836499775519e0; ]; + DD_7(m - 6:m, m - 10:m) = [-0.98986727836499775519e0 0.69173889524389669805e1 -0.20703300533173576296e2 0.34384051724171509761e2 -0.34183525215426919935e2 0.20273644846279336868e2 -0.65319522266785813933e1 0.83355973075426177031e0 0 0 0; 0 -0.97027194845028099974e0 0.67554856358808181957e1 -0.20110258902816885140e2 0.33112901953037622660e2 -0.32392947635031872133e2 0.18405654623600056093e2 -0.55732122985135573041e1 0.77264857229409862775e0 0 0; 0 0 -0.97367953737208690559e0 0.67732979111762465861e1 -0.20116788510903340593e2 0.32902110861795762152e2 -0.31322380032948075240e2 0.18129991367652287552e2 -0.68109221539151178574e1 0.14183700945143243060e1 0; 0 0 0 -0.10718455919447922848e1 0.75637188148594109098e1 -0.22852469704417085459e2 0.37745934806057657407e2 -0.40909341259657284245e2 0.35082323899232605336e2 -0.25816283570514673553e2 0.10257962606384161889e2; 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0; ]; + DD_7 = sparse(DD_7); + + DD_8 = (diag(ones(m - 4, 1), 4) - 8 * diag(ones(m - 3, 1), 3) + 28 * diag(ones(m - 2, 1), 2) - 56 * diag(ones(m - 1, 1), 1) + 70 * diag(ones(m, 1), 0) - 56 * diag(ones(m - 1, 1), -1) + 28 * diag(ones(m - 2, 1), -2) - 8 * diag(ones(m - 3, 1), -3) + diag(ones(m - 4, 1), -4)); + DD_8(1:8, 1:12) = [0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0; 0.11211983054345146839e2 -0.29767555907566203748e2 0.45789440151310044599e2 -0.64530162797168947524e2 0.75491869612115314813e2 -0.60939919211778894557e2 0.30254875259437643639e2 -0.85747647355583382784e1 0.10642345748642342173e1 0 0 0; 0 0.14294303785648149613e1 -0.76427065234692260122e1 0.23888038475126051744e2 -0.50115808052716920383e2 0.65804221723591524304e2 -0.53644769362408908249e2 0.27093191644704986344e2 -0.77894362989766952447e1 0.97783801558437253538e0 0 0; 0 0 0.76035645561041265933e0 -0.63048462739199410204e1 0.24540872831466741457e2 -0.51828716216050995413e2 0.66225803906075245320e2 -0.53627357074178360373e2 0.27021942543523272783e2 -0.77621755876022479980e1 0.97411941507587258409e0 0; 0 0 0 0.82615990808320718845e0 -0.74650882590612358780e1 0.27031526461705782491e2 -0.54693640344683071895e2 0.68768103448343019521e2 -0.55208801421796203457e2 0.27669555809755867922e2 -0.79189382269199820415e1 0.99112262457261614959e0; ]; + DD_8(m - 7:m, m - 11:m) = [0.99112262457261614959e0 -0.79189382269199820415e1 0.27669555809755867922e2 -0.55208801421796203457e2 0.68768103448343019521e2 -0.54693640344683071895e2 0.27031526461705782491e2 -0.74650882590612358780e1 0.82615990808320718845e0 0 0 0; 0 0.97411941507587258409e0 -0.77621755876022479980e1 0.27021942543523272783e2 -0.53627357074178360373e2 0.66225803906075245320e2 -0.51828716216050995413e2 0.24540872831466741457e2 -0.63048462739199410204e1 0.76035645561041265933e0 0 0; 0 0 0.97783801558437253538e0 -0.77894362989766952447e1 0.27093191644704986344e2 -0.53644769362408908249e2 0.65804221723591524304e2 -0.50115808052716920383e2 0.23888038475126051744e2 -0.76427065234692260122e1 0.14294303785648149613e1 0; 0 0 0 0.10642345748642342173e1 -0.85747647355583382784e1 0.30254875259437643639e2 -0.60939919211778894557e2 0.75491869612115314813e2 -0.64530162797168947524e2 0.45789440151310044599e2 -0.29767555907566203748e2 0.11211983054345146839e2; 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0; ]; + DD_8 = sparse(DD_8); + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + %%%% Difference operators %% + D1 = H \ Q; + + % Helper functions for constructing D2(c) + % TODO: Consider changing sparse(diag(...)) to spdiags(....) + + % Minimal 9 point stencil width + function D2 = D2_fun_minimal(c) + % Here we add variable diffusion + C1 = sparse(diag(c)); + C2 = 1/2 * diag(ones(m - 1, 1), -1) + 1/2 * diag(ones(m, 1), 0); C2(1, 2) = 1/2; + C3 = 3/10 * diag(ones(m - 1, 1), -1) + 3/10 * diag(ones(m - 1, 1), 1) + 2/5 * diag(ones(m, 1), 0); C3(1, 3) = 3/10; C3(m, m - 2) = 3/10; + C4 = 1/4 * diag(ones(m - 2, 1), -2) + 1/4 * diag(ones(m - 1, 1), -1) + 1/4 * diag(ones(m - 1, 1), 1) + 1/4 * diag(ones(m, 1), 0); C4(2, 4) = 1/4; C4(1, 3) = 1/4; C4(1, 4) = 1/4; C4(m, m - 2) = 1/4; C4(m, m - 3) = 1/4; + + C2 = sparse(diag(C2 * c)); + C3 = sparse(diag(C3 * c)); + C4 = sparse(diag(C4 * c)); + + % Remainder term added to wide second derivative operator + R = (1/78400 / h) * transpose(DD_8) * C1 * DD_8 + (1/14700 / h) * transpose(DD_7) * C2 * DD_7 + (1/2520 / h) * transpose(DD_6) * C3 * DD_6 + (1/350 / h) * transpose(DD_5) * C4 * DD_5; + + D2 = D1 * C1 * D1 - H \ R; + end + + % Few additional grid point in interior stencil cmp. to minimal + function D2 = D2_fun_nonminimal(c) + % Here we add variable diffusion + C1 = sparse(diag(c)); + C2 = 1/2 * diag(ones(m - 1, 1), -1) + 1/2 * diag(ones(m, 1), 0); C2(1, 2) = 1/2; + C3 = 3/10 * diag(ones(m - 1, 1), -1) + 3/10 * diag(ones(m - 1, 1), 1) + 2/5 * diag(ones(m, 1), 0); C3(1, 3) = 3/10; C3(m, m - 2) = 3/10; + + C2 = sparse(diag(C2 * c)); + C3 = sparse(diag(C3 * c)); + + % Remainder term added to wide second derivative operator + R = (1/78400 / h) * transpose(DD_8) * C1 * DD_8 + (1/14700 / h) * transpose(DD_7) * C2 * DD_7 + (1/2520 / h) * transpose(DD_6) * C3 * DD_6; + D2 = D1 * C1 * D1 - H \ R; + end + + % Wide stencil + function D2 = D2_fun_wide(c) + % Here we add variable diffusion + C1 = sparse(diag(c)); + D2 = D1 * C1 * D1; + end + + switch options.stencil_width + case 'minimal' + D2 = @D2_fun_minimal; + case 'nonminimal' + D2 = @D2_fun_nonminimal; + case 'wide' + D2 = @D2_fun_wide; + otherwise + error('No option %s for stencil width', options.stencil_width) + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%% + + %%%% Artificial dissipation operator %%% + switch options.AD + case 'upwind' + % This is the choice that yield 7th order Upwind + DI = H \ (transpose(DD_4) * DD_4) * (-1/280); + case 'op' + % This choice will preserve the order of the underlying + % Non-dissipative D1 SBP operator + DI = H \ (transpose(DD_5) * DD_5) * (-1 / (5 * 280)); + % Notice that you can use any negative number instead of (-1/(5*280)) + otherwise + error("Artificial dissipation options '%s' not implemented.", option.AD) + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%% +end \ No newline at end of file
--- a/+sbp/D1Nonequidistant.m Fri Mar 13 12:07:00 2020 +0100 +++ b/+sbp/D1Nonequidistant.m Mon Feb 14 11:14:46 2022 +0100 @@ -19,58 +19,53 @@ % 'Accurate' operators are optimized for accuracy % 'Minimal' operators have the smallest possible boundary % closure - - x_l = lim{1}; - x_r = lim{2}; - L = x_r-x_l; - switch option case {'Accurate','accurate','A'} - + [x,h] = sbp.grid.accurateBoundaryOptimizedGrid(lim,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.grid.minimalBoundaryOptimizedGrid(lim,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; obj.e_l = sparse(m,1); obj.e_r = sparse(m,1);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/D2Nonequidistant.m Mon Feb 14 11:14:46 2022 +0100 @@ -0,0 +1,97 @@ +classdef D2Nonequidistant < sbp.OpSet + % Implements the boundary optimized variable coefficient + % second derivative. + % + % The boundary closure uses the first and last rows of the + % boundary-optimized D1 operator, i.e. the operators are + % fully compatible. + properties + H % Norm matrix + HI % H^-1 + D1 % SBP operator approximating first derivative + D2 % SBP operator approximating second derivative + DI % Dissipation operator + e_l % Left boundary operator + e_r % Right boundary operator + d1_l % Left boundary first derivative + d1_r % Right boundary first derivative + m % Number of grid points. + h % Step size + x % grid + borrowing % Struct with borrowing limits for different norm matrices + options % Struct holding options used to create the operator + end + + methods + + function obj = D2Nonequidistant(m, lim, order, options) + % m - number of gridpoints + % lim - cell array holding the limits of the domain + % order - order of the operator + % options - struct holding options used to construct the operator + % struct.stencil_width: {'minimal', 'nonminimal', 'wide'} + % minimal: minimal compatible stencil width (default) + % nonminimal: a few additional stencil points compared to minimal + % wide: wide stencil obtained by applying D1 twice + % struct.AD: {'op', 'upwind'} + % 'op': order-preserving AD (preserving interior stencil order) (default) + % 'upwind': upwind AD (order-1 upwind interior stencil) + % struct.variable_coeffs: {true, false} + % true: obj.D2 is a function handle D2(c) returning a matrix + % for coefficient vector c (default) + % false: obj.D2 is a matrix. + default_arg('options', struct); + default_field(options,'stencil_width','minimal'); + default_field(options,'AD','op'); + default_field(options,'variable_coeffs',true); + [x, h] = sbp.grid.accurateBoundaryOptimizedGrid(lim, m, order); + switch order + case 4 + [obj.H, obj.HI, obj.D1, obj.D2, obj.DI] = sbp.implementations.d2_noneq_variable_4(m, h, options); + case 6 + % [obj.H, obj.HI, obj.D1, D2, M, Q, obj.e_l, obj.e_r, obj.d1_l, obj.d1_r] = sbp.implementations.d2_noneq_6(m,h,2); + % obj.D2 = @(c)D2; + [obj.H, obj.HI, obj.D1, obj.D2, obj.DI] = sbp.implementations.d2_noneq_variable_6(m, h, options); + case 8 + % [obj.H, obj.HI, obj.D1, D2, M, Q, obj.e_l, obj.e_r, obj.d1_l, obj.d1_r] = sbp.implementations.d2_noneq_8(m,h,1); + % obj.D2 = @(c)D2; + [obj.H, obj.HI, obj.D1, obj.D2, obj.DI] = sbp.implementations.d2_noneq_variable_8(m, h, options); + case 10 + [obj.H, obj.HI, obj.D1, obj.D2, obj.DI] = sbp.implementations.d2_noneq_variable_10(m, h, options); + case 12 + [obj.H, obj.HI, obj.D1, obj.D2, obj.DI] = sbp.implementations.d2_noneq_variable_12(m, h, options); + otherwise + error('Invalid operator order %d.', order); + end + + if ~options.variable_coeffs + obj.D2 = obj.D2(ones(m,1)); + end + + % Boundary operators + obj.e_l = sparse(m, 1); obj.e_l(1) = 1; + obj.e_r = sparse(m, 1); obj.e_r(m) = 1; + obj.d1_l = (obj.e_l' * obj.D1)'; + obj.d1_r = (obj.e_r' * obj.D1)'; + + % Borrowing coefficients + obj.borrowing.H11 = obj.H(1, 1) / h; % First element in H/h, + obj.borrowing.M.d1 = obj.H(1, 1) / h; % First element in H/h is borrowing also for M + obj.borrowing.R.delta_D = inf; + + % grid data + obj.x = x; + obj.h = h; + obj.m = m; + + % misc + obj.options = options; + end + + function str = string(obj) + str = [class(obj) '_' num2str(obj.order)]; + end + + end + +end