Mercurial > repos > public > sbplib

comparison +sbp/ordinary12.m @ 255:df3cc9c5dffc operator_remake

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Added ordinary 12th order accurate, with D1*D1 as 2nd derivative.
author Martin Almquist <martin.almquist@it.uu.se>
date Wed, 07 Sep 2016 15:54:41 +0200
parents
children 5714fda0db85
comparison
equal deleted inserted replaced
254:9c97eb0d61bf 255:df3cc9c5dffc
1 function [H, HI, D1, D2, e_1, e_m, M, Q, S_1, S_m] = ordinary12(m,h)
2
3 % D2 = D1*D1, wide!
4
5 H=diag(ones(m,1),0);
6 H(1:15,1:15)=diag([2.880607858916397e-01,...
7 1.560376162339675e+00,...
8 2.403139445479289e-01,...
9 1.531292896754208e+00,...
10 1.245107365551116e+00,...
11 2.403139445479289e-01,...
12 1.218498444843265e+00,...
13 1.730095305997005e+00,...
14 2.403139445479289e-01,...
15 8.526541262207157e-01,...
16 1.839467231356713e+00,...
17 2.403139445479289e-01,...
18 1.350211172437031e+00,...
19 9.139950953929686e-01,...
20 1.008985635023947e+00]);
21
22 H(m-14:m,m-14:m)=rot90(H(1:15,1:15),2);
23
24 H=H*h;
25 HI = inv(H);
26
27 a1 = 6/7; a2 = -15/56; a3 = 5/63; a4 = -1/56; a5 = 1/385; a6 = -1/5544;
28
29 D1=a6*(diag(ones(m-6,1),6)-diag(ones(m-6,1),-6))+...
30 a5*(diag(ones(m-5,1),5)-diag(ones(m-5,1),-5))+...
31 a4*(diag(ones(m-4,1),4)-diag(ones(m-4,1),-4))+...
32 a3*(diag(ones(m-3,1),3)-diag(ones(m-3,1),-3))+...
33 a2*(diag(ones(m-2,1),2)-diag(ones(m-2,1),-2))+...
34 a1*(diag(ones(m-1,1),1)-diag(ones(m-1,1),-1));
35
36
37 D1(1, 1)= -1.735744761135539e+00;
38 D1(1 , 2 )= 2.307154400697361e+00;
39 D1(1, 3)= -8.049541548921420e-02;
40 D1(1 , 4)= -5.054323197649326e-01;
41 D1(1 , 5 )= -2.276940626773821e-01;
42 D1(1 , 6 )= -4.332046513157440e-03;
43 D1(1 , 7 )= 1.922516500571120e-01;
44 D1(1 , 8 )= 1.841426475531248e-01;
45 D1(1 , 9 )= 4.299540500353117e-02;
46 D1(1 , 10 )= -8.425321510779284e-02;
47 D1(1 , 11 )= -1.783432239793461e-01;
48 D1(1 , 12 )= -9.098340284317097e-02;
49 D1(1 , 13 )= 2.989603414709966e-01;
50 D1(1 , 14 )= -1.262666707971657e-01;
51 D1(1 , 15 )= 8.040673525575449e-03;
52 D1(2 , 1 )= -4.259233932680148e-01;
53 D1(2 , 3 )= 1.077970170802666e-01;
54 D1(2 , 4 )= 3.305658398446859e-01;
55 D1(2 , 5 )= 1.023936809322525e-01;
56 D1(2 , 6 )= 1.004615991085861e-02;
57 D1(2 , 7 )= -9.993439452786020e-02;
58 D1(2 , 8 )= -8.071677622127596e-02;
59 D1(2 , 9 )= -2.794842118590027e-02;
60 D1(2 , 10 )= 4.202296437592187e-02;
61 D1(2 , 11 )= 8.084986053360355e-02;
62 D1(2 , 12 )= 4.970793893896007e-02;
63 D1(2 , 13 )= -1.480838886705940e-01;
64 D1(2 , 14 )= 6.449962737470759e-02;
65 D1(2 , 15 )= -5.276215117611654e-03;
66 D1(3 , 1 )= 9.648866897889277e-02;
67 D1(3 , 2 )= -6.999339806925928e-01;
68 D1(3 , 4 )= 6.447106132171313e-01;
69 D1(3 , 5 )= 1.610341017868972e-01;
70 D1(3 , 6 )= -2.926101271261580e-01;
71 D1(3 , 7 )= 2.578035970538076e-02;
72 D1(3 , 8 )= -3.649479605631123e-02;
73 D1(3 , 9 )= 2.039707558852247e-01;
74 D1(3 , 10 )= -2.443704472345096e-02;
75 D1(3 , 11 )= -4.983260374057472e-02;
76 D1(3 , 12 )= -1.982796172639611e-01;
77 D1(3 , 13 )= 2.659121147554237e-01;
78 D1(3 , 14 )= -1.104398155224152e-01;
79 D1(3 , 15 )= 1.413137079651327e-02;
80 D1(4 , 1 )= 9.507993640872410e-02;
81 D1(4 , 2 )= -3.368441515472113e-01;
82 D1(4 , 3 )= -1.011778679849722e-01;
83 D1(4 , 5 )= 1.664700006607916e-01;
84 D1(4 , 6 )= 8.420689674052668e-02;
85 D1(4 , 7 )= 1.234483107384551e-01;
86 D1(4 , 8 )= 3.371285340885326e-02;
87 D1(4 , 9 )= -2.009043511216467e-02;
88 D1(4 , 10 )= -4.311340633236862e-02;
89 D1(4 , 11 )= -3.824719371220764e-02;
90 D1(4 , 12 )= -2.775022659375772e-03;
91 D1(4 , 13 )= 8.248678009703486e-02;
92 D1(4 , 14 )= -5.169933586327551e-02;
93 D1(4 , 15 )= 8.542635157189999e-03;
94 D1(5 , 1 )= 5.267797175762052e-02;
95 D1(5 , 2 )= -1.283203869171410e-01;
96 D1(5 , 3 )= -3.108064515385221e-02;
97 D1(5 , 4 )= -2.047328098663259e-01;
98 D1(5 , 6 )= 1.528377576445926e-02;
99 D1(5 , 7 )= 1.842140733907868e-01;
100 D1(5 , 8 )= 1.776491663357457e-01;
101 D1(5 , 9 )= 4.220032524144508e-02;
102 D1(5 , 10 )= -3.611067700750285e-02;
103 D1(5 , 11 )= -1.090730089234403e-01;
104 D1(5 , 12 )= -5.398262792345335e-02;
105 D1(5 , 13 )= 1.325672256468355e-01;
106 D1(5 , 14 )= -3.865175590539601e-02;
107 D1(5 , 15 )= -2.640626439781275e-03;
108 D1(6 , 1 )= 5.192760351242898e-03;
109 D1(6 , 2 )= -6.523045708997477e-02;
110 D1(6 , 3 )= 2.926101271261579e-01;
111 D1(6 , 4 )= -5.365707057867646e-01;
112 D1(6 , 5 )= -7.918783828195401e-02;
113 D1(6 , 7 )= 1.455199260610844e-01;
114 D1(6 , 8 )= 4.392627009892030e-01;
115 D1(6 , 9 )= -2.877704902282471e-01;
116 D1(6 , 10 )= 6.246304940524665e-02;
117 D1(6 , 11 )= -8.757523368284212e-03;
118 D1(6 , 12 )= 2.373907173856874e-01;
119 D1(6 , 13 )= -3.887671913712637e-01;
120 D1(6 , 14 )= 2.338336212010801e-01;
121 D1(6 , 15 )= -4.998869639321427e-02;
122 D1( 7 , 1 )= -4.544951340626429e-02;
123 D1( 7 , 2 )= 1.279732835762292e-01;
124 D1( 7 , 3 )= -5.084438112238577e-03;
125 D1( 7 , 4 )= -1.551380899582667e-01;
126 D1( 7 , 5 )= -1.882368423100825e-01;
127 D1( 7 , 6 )= -2.869964060279157e-02;
128 D1( 7 , 8 )= 1.501531969717210e-01;
129 D1( 7 , 9 )= 4.398683557974375e-02;
130 D1( 7 , 10 )= 9.034991722085359e-02;
131 D1( 7 , 11 )= 6.874263504449905e-02;
132 D1( 7 , 12 )= 1.191424140454980e-02;
133 D1( 7 , 13 )= -1.387149470006835e-01;
134 D1( 7, 14 )= 8.081158783540948e-02;
135 D1( 7 , 15 )= -1.260822624267884e-02;
136 D1( 8 , 1 )= -3.065974203065793e-02;
137 D1( 8 , 2 )= 7.279861004189266e-02;
138 D1( 8 , 3 )= 5.069205358435631e-03;
139 D1( 8 , 4 )= -2.983896481040582e-02;
140 D1( 8 , 5 )= -1.278497691554543e-01;
141 D1( 8 , 6 )= -6.101453024095747e-02;
142 D1( 8 , 7 )= -1.057522301598588e-01;
143 D1( 8 , 9 )= 3.518098475279482e-02;
144 D1( 8 , 10 )= 1.267951379582364e-01;
145 D1( 8 , 11 )= 1.709119775233402e-01;
146 D1( 8 , 12 )= 4.871305182138817e-02;
147 D1( 8 , 13 )= -1.381315738947210e-01;
148 D1( 8 , 14 )= 2.531339542383916e-02;
149 D1( 8 , 15 )= 8.464447412128264e-03;
150 D1( 9 , 1 )= -5.153795872455654e-02;
151 D1( 9 , 2 )= 1.814711596347261e-01;
152 D1( 9 , 3 )= -2.039707558852247e-01;
153 D1( 9 , 4 )= 1.280172927036423e-01;
154 D1( 9 , 5 )= -2.186470530689342e-01;
155 D1( 9 , 6 )= 2.877704902282471e-01;
156 D1( 9 , 7 )= -2.230327950727989e-01;
157 D1( 9 , 8 )= -2.532789210200124e-01;
158 D1( 9 , 10 )= -2.487018120584765e-02;
159 D1( 9 , 11 )= 5.656432879788119e-01;
160 D1( 9 , 12 )= -3.179024054582303e-01;
161 D1( 9 , 13 )= 3.080272260545759e-01;
162 D1( 9 , 14 )= -2.449533574878925e-01;
163 D1( 9 , 15 )= 6.726397132349360e-02;
164 D1(10 , 1 )= 2.846411764336633e-02;
165 D1(10 , 2 )= -7.690296670899374e-02;
166 D1(10 , 3 )= 6.887391299700927e-03;
167 D1(10 , 4 )= 7.742794040563218e-02;
168 D1(10 , 5 )= 5.273142829480713e-02;
169 D1(10 , 6 )= -1.760472544430208e-02;
170 D1(10 , 7 )= -1.291159336943500e-01;
171 D1(10 , 8 )= -2.572762697778848e-01;
172 D1(10 , 9 )= 7.009467453925059e-03;
173 D1(10 , 11 )= 1.983435406810274e-01;
174 D1(10 , 12 )= -3.844771590180474e-02;
175 D1(10 , 13 )= 2.943472385247488e-01;
176 D1(10 , 14 )= -1.771193581329739e-01;
177 D1(10 , 15 )= 3.146739090039766e-02;
178 D1(10 , 16 )= -2.115455432962849e-04;
179 D1(11 , 1 )= 2.792856995883971e-02;
180 D1(11 , 2 )= -6.858300759838749e-02;
181 D1(11 , 3 )= 6.510292419375574e-03;
182 D1(11 , 4 )= 3.183946691406334e-02;
183 D1(11 , 5 )= 7.382985925399284e-02;
184 D1(11 , 6 )= 1.144111158506918e-03;
185 D1(11 , 7 )= -4.553644254612264e-02;
186 D1(11 , 8 )= -1.607498111470597e-01;
187 D1(11 , 9 )= -7.389746738841917e-02;
188 D1(11 , 10 )= -9.193881548309495e-02;
189 D1(11 , 12 )= 4.858403641270599e-02;
190 D1(11 , 13 )= 2.953396945760745e-01;
191 D1(11 , 14 )= -1.897059300798801e-02;
192 D1(11 , 15 )= -2.681387584075025e-02;
193 D1(11 , 16 )= 1.412040700223218e-03;
194 D1(11 , 17 )= -9.805838195994568e-05;
195 D1(12 , 1 )= 1.090604649488924e-01;
196 D1(12 , 2 )= -3.227573129195589e-01;
197 D1(12 , 3 )= 1.982796172639611e-01;
198 D1(12 , 4 )= 1.768258806049677e-02;
199 D1(12 , 5 )= 2.796931645632894e-01;
200 D1(12 , 6 )= -2.373907173856874e-01;
201 D1(12 , 7 )= -6.041049615427440e-02;
202 D1(12 , 8 )= -3.507005074362792e-01;
203 D1(12 , 9 )= 3.179024054582303e-01;
204 D1(12 , 10 )= 1.364157359619944e-01;
205 D1(12 , 11 )= -3.718833008893085e-01;
206 D1(12 , 13 )= -1.473587211078091e-01;
207 D1(12 , 14 )= 5.107219886498376e-01;
208 D1(12 , 15 )= -1.500513994956780e-02;
209 D1(12 , 16 )= -7.430756001585824e-02;
210 D1(12 , 17 )= 1.080837236594302e-02;
211 D1(12 , 18 )= -7.505814143015984e-04;
212 D1(13 , 1 )= -6.378169035524298e-02;
213 D1(13 , 2 )= 1.711336527390002e-01;
214 D1(13 , 3 )= -4.732769992164881e-02;
215 D1(13 , 4 )= -9.354938176872968e-02;
216 D1(13 , 5 )= -1.222478620034153e-01;
217 D1(13 , 6 )= 6.919375218960795e-02;
218 D1(13 , 7 )= 1.251833421669687e-01;
219 D1(13 , 8 )= 1.769951193441043e-01;
220 D1(13 , 9 )= -5.482345223653074e-02;
221 D1(13 , 10 )= -1.858793591648383e-01;
222 D1(13 , 11 )= -4.023575729350775e-01;
223 D1(13 , 12 )= 2.622727189335804e-02;
224 D1(13 , 14 )= 4.511356466215911e-01;
225 D1(13 , 15 )= -9.724618589400275e-02;
226 D1(13 , 16 )= 5.877975311212341e-02;
227 D1(13 , 17 )= -1.322544445022777e-02;
228 D1(13 , 18 )= 1.923701010942221e-03;
229 D1(13 , 19 )= -1.335903479820987e-04;
230 D1(14 , 1 )= 3.979504551511218e-02;
231 D1(14 , 2 )= -1.101140274631495e-01;
232 D1(14 , 3 )= 2.903760407152547e-02;
233 D1(14 , 4 )= 8.661624791356937e-02;
234 D1(14 , 5 )= 5.265409651744463e-02;
235 D1(14 , 6 )= -6.148116129069339e-02;
236 D1(14 , 7 )= -1.077344885099467e-01;
237 D1(14 , 8 )= -4.791555974685060e-02;
238 D1(14 , 9 )= 6.440483965930399e-02;
239 D1(14 , 10 )= 1.652323434850750e-01;
240 D1(14 , 11 )= 3.817939983867793e-02;
241 D1(14 , 12 )= -1.342825757801646e-01;
242 D1(14 , 13 )= -6.664460164211102e-01;
243 D1(14 , 15 )= 8.751759953804069e-01;
244 D1(14 , 16 )= -2.930619039503475e-01;
245 D1(14 , 17 )= 8.683315672602888e-02;
246 D1(14 , 18 )= -1.953746026335650e-02;
247 D1(14 , 19 )= 2.841812401942763e-03;
248 D1(14 , 20 )= -1.973480834682474e-04;
249 D1(15 , 1 )= -2.295575530984041e-03;
250 D1(15 , 2 )= 8.159561455701054e-03;
251 D1(15 , 3 )= -3.365722305748102e-03;
252 D1(15 , 4 )= -1.296477975670818e-02;
253 D1(15 , 5 )= 3.258582992375954e-03;
254 D1(15 , 6 )= 1.190599786168119e-02;
255 D1(15 , 7 )= 1.522628621820923e-02;
256 D1(15 , 8 )= -1.451388426876270e-02;
257 D1(15 , 9 )= -1.602051576712872e-02;
258 D1(15 , 10 )= -2.659185598017704e-02;
259 D1(15 , 11 )= 4.888399224192811e-02;
260 D1(15 , 12 )= 3.573831226733737e-03;
261 D1(15 , 13 )= 1.301335540499094e-01;
262 D1(15 , 14 )= -7.927829095053140e-01;
263 D1(15 , 16 )= 8.495094750506673e-01;
264 D1(15 , 17 )= -2.654717109533336e-01;
265 D1(15 , 18 )= 7.865828472691364e-02;
266 D1(15 , 19 )= -1.769811406355557e-02;
267 D1(15 , 20 )= 2.574271136517173e-03;
268 D1(15 , 21 )= -1.787688289248037e-04;
269
270
271 D1(m-14:m,m-20:m)=rot90( -D1(1:15,1:21),2);
272 D1=D1/h;
273
274
275 % Wide 2nd derivative
276
277 D2=D1*D1;
278
279 e_1 = zeros(m,1);
280 e_1(1)= 1;
281 e_m = zeros(m,1);
282 e_m(end)= 1;
283
284 S_1 = (e_1'*D1)';
285 S_m = (e_m'*D1);
286
287 Q = H*D1-(-e_1*e_1' + e_m*e_m');
288 M = -(H*D2-(-e_1*S_1' + e_m*S_m'));
289 end