Mercurial > repos > public > sbplib

comparison +sbp/+implementations/d1_noneq_12.m @ 261:6009f2712d13 operator_remake

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