Mercurial > repos > public > sbplib
comparison +sbp/+implementations/d1_noneq_6.m @ 1286:4cb627c7fb90 feature/boundary_optimized_grids
Make D1Nonequidistant use the grid generation functions accurate/minimalBoundaryOptimizedGrid and remove grid generation from +implementations
author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
---|---|
date | Wed, 01 Jul 2020 13:43:32 +0200 |
parents | f7ac3cd6eeaa |
children |
comparison
equal
deleted
inserted
replaced
1285:6b68f939d023 | 1286:4cb627c7fb90 |
---|---|
1 function [D1,H,x,h] = d1_noneq_6(N,L) | 1 function [D1,H] = d1_noneq_6(N,h) |
2 | 2 |
3 % L: Domain length | |
4 % N: Number of grid points | 3 % N: Number of grid points |
5 if(nargin < 2) | |
6 L = 1; | |
7 end | |
8 | |
9 if(N<12) | 4 if(N<12) |
10 error('Operator requires at least 12 grid points'); | 5 error('Operator requires at least 12 grid points'); |
11 end | 6 end |
12 | 7 |
13 % BP: Number of boundary points | 8 % BP: Number of boundary points |
14 % m: Number of nonequidistant spacings | |
15 % order: Accuracy of interior stencil | |
16 BP = 6; | 9 BP = 6; |
17 m = 3; | |
18 order = 6; | |
19 | |
20 %%%% Non-equidistant grid points %%%%% | |
21 x0 = 0.0000000000000e+00; | |
22 x1 = 4.4090263368623e-01; | |
23 x2 = 1.2855984345073e+00; | |
24 x3 = 2.2638953951239e+00; | |
25 x4 = 3.2638953951239e+00; | |
26 x5 = 4.2638953951239e+00; | |
27 x6 = 5.2638953951239e+00; | |
28 | |
29 xb = sparse(m+1,1); | |
30 for i = 0:m | |
31 xb(i+1) = eval(['x' num2str(i)]); | |
32 end | |
33 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
34 | |
35 %%%% Compute h %%%%%%%%%% | |
36 h = L/(2*xb(end) + N-1-2*m); | |
37 %%%%%%%%%%%%%%%%%%%%%%%%% | |
38 | |
39 %%%% Define grid %%%%%%%% | |
40 x = h*[xb; linspace(xb(end)+1,L/h-xb(end)-1,N-2*(m+1))'; L/h-flip(xb) ]; | |
41 %%%%%%%%%%%%%%%%%%%%%%%%% | |
42 | 10 |
43 %%%% Norm matrix %%%%%%%% | 11 %%%% Norm matrix %%%%%%%% |
44 P = sparse(BP,1); | 12 P = sparse(BP,1); |
45 %#ok<*NASGU> | 13 %#ok<*NASGU> |
46 P0 = 1.3030223027124e-01; | 14 P0 = 1.3030223027124e-01; |
59 H(end-BP+1:end) = flip(P); | 27 H(end-BP+1:end) = flip(P); |
60 H = spdiags(h*H,0,N,N); | 28 H = spdiags(h*H,0,N,N); |
61 %%%%%%%%%%%%%%%%%%%%%%%%% | 29 %%%%%%%%%%%%%%%%%%%%%%%%% |
62 | 30 |
63 %%%% Q matrix %%%%%%%%%%% | 31 %%%% Q matrix %%%%%%%%%%% |
64 | |
65 % interior stencil | 32 % interior stencil |
66 switch order | 33 order = 6; |
67 case 2 | 34 d = [-1/60,3/20,-3/4,0,3/4,-3/20,1/60]; |
68 d = [-1/2,0,1/2]; | |
69 case 4 | |
70 d = [1/12,-2/3,0,2/3,-1/12]; | |
71 case 6 | |
72 d = [-1/60,3/20,-3/4,0,3/4,-3/20,1/60]; | |
73 case 8 | |
74 d = [1/280,-4/105,1/5,-4/5,0,4/5,-1/5,4/105,-1/280]; | |
75 case 10 | |
76 d = [-1/1260,5/504,-5/84,5/21,-5/6,0,5/6,-5/21,5/84,-5/504,1/1260]; | |
77 case 12 | |
78 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]; | |
79 end | |
80 d = repmat(d,N,1); | 35 d = repmat(d,N,1); |
81 Q = spdiags(d,-order/2:order/2,N,N); | 36 Q = spdiags(d,-order/2:order/2,N,N); |
82 | 37 |
83 % Boundaries | 38 % Boundaries |
84 Q0_0 = -5.0000000000000e-01; | 39 Q0_0 = -5.0000000000000e-01; |