Mercurial > repos > public > sbplib

comparison +sbp/+implementations/d1_noneq_12.m @ 423:a2cb0d4f4a02 feature/grids

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge in default.
author Jonatan Werpers <jonatan@werpers.com>
date Tue, 07 Feb 2017 15:47:51 +0100
parents f7ac3cd6eeaa
children 4cb627c7fb90
comparison
equal deleted inserted replaced
218:da058ce66876 423:a2cb0d4f4a02
1 function [D1,H,x,h] = d1_noneq_12(N,L)
2
3 % L: Domain length
4 % N: Number of grid points
5 if(nargin < 2)
6 L = 1;
7 end
8
9 if(N<24)
10 error('Operator requires at least 24 grid points');
11 end
12
13 % BP: Number of boundary points
14 % m: Number of nonequidistant spacings
15 % order: Accuracy of interior stencil
16 BP = 12;
17 m = 6;
18 order = 12;
19
20 %%%% Non-equidistant grid points %%%%%
21 x0 = 0.0000000000000e+00;
22 x1 = 3.6098032343909e-01;
23 x2 = 1.1634317168086e+00;
24 x3 = 2.2975905356987e+00;
25 x4 = 3.6057529790929e+00;
26 x5 = 4.8918275675510e+00;
27 x6 = 6.0000000000000e+00;
28 x7 = 7.0000000000000e+00;
29 x8 = 8.0000000000000e+00;
30 x9 = 9.0000000000000e+00;
31 x10 = 1.0000000000000e+01;
32 x11 = 1.1000000000000e+01;
33 x12 = 1.2000000000000e+01;
34
35 xb = sparse(m+1,1);
36 for i = 0:m
37 xb(i+1) = eval(['x' num2str(i)]);
38 end
39 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
40
41 %%%% Compute h %%%%%%%%%%
42 h = L/(2*xb(end) + N-1-2*m);
43 %%%%%%%%%%%%%%%%%%%%%%%%%
44
45 %%%% Define grid %%%%%%%%
46 x = h*[xb; linspace(xb(end)+1,L/h-xb(end)-1,N-2*(m+1))'; L/h-flip(xb) ];
47 %%%%%%%%%%%%%%%%%%%%%%%%%
48
49 %%%% Norm matrix %%%%%%%%
50 P = sparse(BP,1);
51 %#ok<*NASGU>
52 P0 = 1.0000000000011e-01;
53 P1 = 5.9616216757547e-01;
54 P2 = 9.9065699844442e-01;
55 P3 = 1.2512548713913e+00;
56 P4 = 1.3316678989403e+00;
57 P5 = 1.2093375037721e+00;
58 P6 = 1.0236491595704e+00;
59 P7 = 9.9685258909811e-01;
60 P8 = 1.0004766563923e+00;
61 P9 = 9.9993617879146e-01;
62 P10 = 1.0000063122914e+00;
63 P11 = 9.9999966373260e-01;
64
65 for i = 0:BP-1
66 P(i+1) = eval(['P' num2str(i)]);
67 end
68
69 H = ones(N,1);
70 H(1:BP) = P;
71 H(end-BP+1:end) = flip(P);
72 H = spdiags(h*H,0,N,N);
73 %%%%%%%%%%%%%%%%%%%%%%%%%
74
75 %%%% Q matrix %%%%%%%%%%%
76
77 % interior stencil
78 switch order
79 case 2
80 d = [-1/2,0,1/2];
81 case 4
82 d = [1/12,-2/3,0,2/3,-1/12];
83 case 6
84 d = [-1/60,3/20,-3/4,0,3/4,-3/20,1/60];
85 case 8
86 d = [1/280,-4/105,1/5,-4/5,0,4/5,-1/5,4/105,-1/280];
87 case 10
88 d = [-1/1260,5/504,-5/84,5/21,-5/6,0,5/6,-5/21,5/84,-5/504,1/1260];
89 case 12
90 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];
91 end
92 d = repmat(d,N,1);
93 Q = spdiags(d,-order/2:order/2,N,N);
94
95 % Boundaries
96 Q0_0 = -5.0000000000000e-01;
97 Q0_1 = 6.7597560728423e-01;
98 Q0_2 = -2.6859785384416e-01;
99 Q0_3 = 1.4850302678903e-01;
100 Q0_4 = -8.7976689586154e-02;
101 Q0_5 = 4.1833336322613e-02;
102 Q0_6 = -2.2216684976993e-03;
103 Q0_7 = -1.5910034062022e-02;
104 Q0_8 = 1.1296706376589e-02;
105 Q0_9 = -3.1823678285130e-03;
106 Q0_10 = 2.4843594063649e-04;
107 Q0_11 = 3.1501105449828e-05;
108 Q0_12 = 0.0000000000000e+00;
109 Q0_13 = 0.0000000000000e+00;
110 Q0_14 = 0.0000000000000e+00;
111 Q0_15 = 0.0000000000000e+00;
112 Q0_16 = 0.0000000000000e+00;
113 Q0_17 = 0.0000000000000e+00;
114 Q1_0 = -6.7597560728423e-01;
115 Q1_1 = 0.0000000000000e+00;
116 Q1_2 = 9.5424013647146e-01;
117 Q1_3 = -4.3389334603464e-01;
118 Q1_4 = 2.4285669347653e-01;
119 Q1_5 = -1.1443465137214e-01;
120 Q1_6 = 8.5942765682435e-03;
121 Q1_7 = 4.0290424215772e-02;
122 Q1_8 = -2.9396383714543e-02;
123 Q1_9 = 8.5601827834256e-03;
124 Q1_10 = -7.8128092862319e-04;
125 Q1_11 = -6.0444181254875e-05;
126 Q1_12 = 0.0000000000000e+00;
127 Q1_13 = 0.0000000000000e+00;
128 Q1_14 = 0.0000000000000e+00;
129 Q1_15 = 0.0000000000000e+00;
130 Q1_16 = 0.0000000000000e+00;
131 Q1_17 = 0.0000000000000e+00;
132 Q2_0 = 2.6859785384416e-01;
133 Q2_1 = -9.5424013647146e-01;
134 Q2_2 = 0.0000000000000e+00;
135 Q2_3 = 9.7065114311923e-01;
136 Q2_4 = -4.3205328628292e-01;
137 Q2_5 = 1.9549970932735e-01;
138 Q2_6 = -2.4406885385172e-02;
139 Q2_7 = -5.5737279079895e-02;
140 Q2_8 = 4.3772338637753e-02;
141 Q2_9 = -1.3727655130726e-02;
142 Q2_10 = 1.6271304373071e-03;
143 Q2_11 = 1.7066984372933e-05;
144 Q2_12 = 0.0000000000000e+00;
145 Q2_13 = 0.0000000000000e+00;
146 Q2_14 = 0.0000000000000e+00;
147 Q2_15 = 0.0000000000000e+00;
148 Q2_16 = 0.0000000000000e+00;
149 Q2_17 = 0.0000000000000e+00;
150 Q3_0 = -1.4850302678903e-01;
151 Q3_1 = 4.3389334603464e-01;
152 Q3_2 = -9.7065114311923e-01;
153 Q3_3 = 0.0000000000000e+00;
154 Q3_4 = 9.5375878629204e-01;
155 Q3_5 = -3.6113954384951e-01;
156 Q3_6 = 6.9749289223875e-02;
157 Q3_7 = 6.5161366516465e-02;
158 Q3_8 = -6.0325702283960e-02;
159 Q3_9 = 2.1188913621662e-02;
160 Q3_10 = -3.2632650250470e-03;
161 Q3_11 = 1.3097937809499e-04;
162 Q3_12 = 0.0000000000000e+00;
163 Q3_13 = 0.0000000000000e+00;
164 Q3_14 = 0.0000000000000e+00;
165 Q3_15 = 0.0000000000000e+00;
166 Q3_16 = 0.0000000000000e+00;
167 Q3_17 = 0.0000000000000e+00;
168 Q4_0 = 8.7976689586154e-02;
169 Q4_1 = -2.4285669347653e-01;
170 Q4_2 = 4.3205328628292e-01;
171 Q4_3 = -9.5375878629204e-01;
172 Q4_4 = 0.0000000000000e+00;
173 Q4_5 = 8.8676146394834e-01;
174 Q4_6 = -2.1292503103800e-01;
175 Q4_7 = -4.6037018833218e-02;
176 Q4_8 = 7.4338719466734e-02;
177 Q4_9 = -3.1217656663809e-02;
178 Q4_10 = 6.1239492854797e-03;
179 Q4_11 = -4.5892226603067e-04;
180 Q4_12 = 0.0000000000000e+00;
181 Q4_13 = 0.0000000000000e+00;
182 Q4_14 = 0.0000000000000e+00;
183 Q4_15 = 0.0000000000000e+00;
184 Q4_16 = 0.0000000000000e+00;
185 Q4_17 = 0.0000000000000e+00;
186 Q5_0 = -4.1833336322613e-02;
187 Q5_1 = 1.1443465137214e-01;
188 Q5_2 = -1.9549970932735e-01;
189 Q5_3 = 3.6113954384951e-01;
190 Q5_4 = -8.8676146394834e-01;
191 Q5_5 = 0.0000000000000e+00;
192 Q5_6 = 7.7461223007026e-01;
193 Q5_7 = -1.0609547334165e-01;
194 Q5_8 = -4.4853791547749e-02;
195 Q5_9 = 3.2436468405486e-02;
196 Q5_10 = -8.4387621360184e-03;
197 Q5_11 = 8.5964292632428e-04;
198 Q5_12 = 0.0000000000000e+00;
199 Q5_13 = 0.0000000000000e+00;
200 Q5_14 = 0.0000000000000e+00;
201 Q5_15 = 0.0000000000000e+00;
202 Q5_16 = 0.0000000000000e+00;
203 Q5_17 = 0.0000000000000e+00;
204 Q6_0 = 2.2216684976993e-03;
205 Q6_1 = -8.5942765682435e-03;
206 Q6_2 = 2.4406885385172e-02;
207 Q6_3 = -6.9749289223875e-02;
208 Q6_4 = 2.1292503103800e-01;
209 Q6_5 = -7.7461223007026e-01;
210 Q6_6 = 0.0000000000000e+00;
211 Q6_7 = 7.4758103262966e-01;
212 Q6_8 = -1.5730779067906e-01;
213 Q6_9 = 2.6517620342970e-02;
214 Q6_10 = -4.3175367549700e-03;
215 Q6_11 = 1.1092605832824e-03;
216 Q6_12 = -1.8037518037522e-04;
217 Q6_13 = 0.0000000000000e+00;
218 Q6_14 = 0.0000000000000e+00;
219 Q6_15 = 0.0000000000000e+00;
220 Q6_16 = 0.0000000000000e+00;
221 Q6_17 = 0.0000000000000e+00;
222 Q7_0 = 1.5910034062022e-02;
223 Q7_1 = -4.0290424215772e-02;
224 Q7_2 = 5.5737279079895e-02;
225 Q7_3 = -6.5161366516465e-02;
226 Q7_4 = 4.6037018833218e-02;
227 Q7_5 = 1.0609547334165e-01;
228 Q7_6 = -7.4758103262966e-01;
229 Q7_7 = 0.0000000000000e+00;
230 Q7_8 = 8.0975719267918e-01;
231 Q7_9 = -2.3568822398349e-01;
232 Q7_10 = 6.9373143801571e-02;
233 Q7_11 = -1.6606121869177e-02;
234 Q7_12 = 2.5974025974031e-03;
235 Q7_13 = -1.8037518037522e-04;
236 Q7_14 = 0.0000000000000e+00;
237 Q7_15 = 0.0000000000000e+00;
238 Q7_16 = 0.0000000000000e+00;
239 Q7_17 = 0.0000000000000e+00;
240 Q8_0 = -1.1296706376589e-02;
241 Q8_1 = 2.9396383714543e-02;
242 Q8_2 = -4.3772338637753e-02;
243 Q8_3 = 6.0325702283960e-02;
244 Q8_4 = -7.4338719466734e-02;
245 Q8_5 = 4.4853791547749e-02;
246 Q8_6 = 1.5730779067906e-01;
247 Q8_7 = -8.0975719267918e-01;
248 Q8_8 = 0.0000000000000e+00;
249 Q8_9 = 8.4765775072084e-01;
250 Q8_10 = -2.6369594097148e-01;
251 Q8_11 = 7.8759594625702e-02;
252 Q8_12 = -1.7857142857146e-02;
253 Q8_13 = 2.5974025974031e-03;
254 Q8_14 = -1.8037518037522e-04;
255 Q8_15 = 0.0000000000000e+00;
256 Q8_16 = 0.0000000000000e+00;
257 Q8_17 = 0.0000000000000e+00;
258 Q9_0 = 3.1823678285130e-03;
259 Q9_1 = -8.5601827834256e-03;
260 Q9_2 = 1.3727655130726e-02;
261 Q9_3 = -2.1188913621662e-02;
262 Q9_4 = 3.1217656663809e-02;
263 Q9_5 = -3.2436468405486e-02;
264 Q9_6 = -2.6517620342970e-02;
265 Q9_7 = 2.3568822398349e-01;
266 Q9_8 = -8.4765775072084e-01;
267 Q9_9 = 0.0000000000000e+00;
268 Q9_10 = 8.5631774953989e-01;
269 Q9_11 = -2.6769768119702e-01;
270 Q9_12 = 7.9365079365093e-02;
271 Q9_13 = -1.7857142857146e-02;
272 Q9_14 = 2.5974025974031e-03;
273 Q9_15 = -1.8037518037522e-04;
274 Q9_16 = 0.0000000000000e+00;
275 Q9_17 = 0.0000000000000e+00;
276 Q10_0 = -2.4843594063649e-04;
277 Q10_0 = -2.4843594063649e-04;
278 Q10_1 = 7.8128092862319e-04;
279 Q10_1 = 7.8128092862319e-04;
280 Q10_2 = -1.6271304373071e-03;
281 Q10_2 = -1.6271304373071e-03;
282 Q10_3 = 3.2632650250470e-03;
283 Q10_3 = 3.2632650250470e-03;
284 Q10_4 = -6.1239492854797e-03;
285 Q10_4 = -6.1239492854797e-03;
286 Q10_5 = 8.4387621360184e-03;
287 Q10_5 = 8.4387621360184e-03;
288 Q10_6 = 4.3175367549700e-03;
289 Q10_6 = 4.3175367549700e-03;
290 Q10_7 = -6.9373143801571e-02;
291 Q10_7 = -6.9373143801571e-02;
292 Q10_8 = 2.6369594097148e-01;
293 Q10_8 = 2.6369594097148e-01;
294 Q10_9 = -8.5631774953989e-01;
295 Q10_9 = -8.5631774953989e-01;
296 Q10_10 = 0.0000000000000e+00;
297 Q10_10 = 0.0000000000000e+00;
298 Q10_11 = 8.5712580212095e-01;
299 Q10_11 = 8.5712580212095e-01;
300 Q10_12 = -2.6785714285718e-01;
301 Q10_12 = -2.6785714285718e-01;
302 Q10_13 = 7.9365079365093e-02;
303 Q10_13 = 7.9365079365093e-02;
304 Q10_14 = -1.7857142857146e-02;
305 Q10_14 = -1.7857142857146e-02;
306 Q10_15 = 2.5974025974031e-03;
307 Q10_15 = 2.5974025974031e-03;
308 Q10_16 = -1.8037518037522e-04;
309 Q10_16 = -1.8037518037522e-04;
310 Q10_17 = 0.0000000000000e+00;
311 Q10_17 = 0.0000000000000e+00;
312 Q11_0 = -3.1501105449828e-05;
313 Q11_0 = -3.1501105449828e-05;
314 Q11_1 = 6.0444181254875e-05;
315 Q11_1 = 6.0444181254875e-05;
316 Q11_2 = -1.7066984372933e-05;
317 Q11_2 = -1.7066984372933e-05;
318 Q11_3 = -1.3097937809499e-04;
319 Q11_3 = -1.3097937809499e-04;
320 Q11_4 = 4.5892226603067e-04;
321 Q11_4 = 4.5892226603067e-04;
322 Q11_5 = -8.5964292632428e-04;
323 Q11_5 = -8.5964292632428e-04;
324 Q11_6 = -1.1092605832824e-03;
325 Q11_6 = -1.1092605832824e-03;
326 Q11_7 = 1.6606121869177e-02;
327 Q11_7 = 1.6606121869177e-02;
328 Q11_8 = -7.8759594625702e-02;
329 Q11_8 = -7.8759594625702e-02;
330 Q11_9 = 2.6769768119702e-01;
331 Q11_9 = 2.6769768119702e-01;
332 Q11_10 = -8.5712580212095e-01;
333 Q11_10 = -8.5712580212095e-01;
334 Q11_11 = 0.0000000000000e+00;
335 Q11_11 = 0.0000000000000e+00;
336 Q11_12 = 8.5714285714289e-01;
337 Q11_12 = 8.5714285714289e-01;
338 Q11_13 = -2.6785714285718e-01;
339 Q11_13 = -2.6785714285718e-01;
340 Q11_14 = 7.9365079365093e-02;
341 Q11_14 = 7.9365079365093e-02;
342 Q11_15 = -1.7857142857146e-02;
343 Q11_15 = -1.7857142857146e-02;
344 Q11_16 = 2.5974025974031e-03;
345 Q11_16 = 2.5974025974031e-03;
346 Q11_17 = -1.8037518037522e-04;
347 Q11_17 = -1.8037518037522e-04;
348 for i = 1:BP
349 for j = 1:BP
350 Q(i,j) = eval(['Q' num2str(i-1) '_' num2str(j-1)]);
351 Q(N+1-i,N+1-j) = -eval(['Q' num2str(i-1) '_' num2str(j-1)]);
352 end
353 end
354 %%%%%%%%%%%%%%%%%%%%%%%%%%%
355
356 %%%% Difference operator %%
357 D1 = H\Q;
358 %%%%%%%%%%%%%%%%%%%%%%%%%%%