1 /*
2 * This source code is part of
3 *
4 * E R K A L E
5 * -
6 * DFT from Hel
7 *
8 * Written by Susi Lehtola, 2011
9 * Copyright (c) 2011, Susi Lehtola
10 *
11 * This program is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU General Public License
13 * as published by the Free Software Foundation; either version 2
14 * of the License, or (at your option) any later version.
15 *
16 * This file contains the Lobatto quadrature routine ported from quadrule.
17 * The original license was LGPL, this version is relicensed to GPLv2+.
18 */
19
20
21 #include <cfloat>
22 #include <cmath>
23 #include <cstdio>
24 #include <sstream>
25 #include <stdexcept>
26
27 #include "lobatto.h"
28 #include "mathf.h"
29
lobatto_set(int order,std::vector<double> & xtab,std::vector<double> & weight)30 void lobatto_set(int order, std::vector<double> & xtab, std::vector<double> & weight)
31
32 /******************************************************************************/
33 /*
34 Purpose:
35
36 LOBATTO_SET sets abscissas and weights for Lobatto quadrature.
37
38 Discussion:
39
40 The integral:
41
42 Integral ( -1 <= X <= 1 ) F(X) dX
43
44 The quadrature rule:
45
46 Sum ( 1 <= I <= ORDER ) WEIGHt[I) * F ( XTAb[I) )
47
48 The quadrature rule will integrate exactly all polynomials up to
49 X**(2*ORDER-3).
50
51 The Lobatto rule is distinguished by the fact that both endpoints
52 (-1 and 1) are always abscissas.
53
54 Licensing:
55
56 This code is distributed under the GNU LGPL license.
57
58 Modified:
59
60 30 April 2006
61
62 Author:
63
64 John Burkardt
65
66 Reference:
67
68 Milton Abramowitz, Irene Stegun,
69 Handbook of Mathematical Functions,
70 National Bureau of Standards, 1964,
71 ISBN: 0-486-61272-4,
72 LC: QA47.A34.
73
74 Arthur Stroud, Don Secrest,
75 Gaussian Quadrature Formulas,
76 Prentice Hall, 1966,
77 LC: QA299.4G3S7.
78
79 Daniel Zwillinger, editor,
80 CRC Standard Mathematical Tables and Formulae,
81 30th Edition,
82 CRC Press, 1996,
83 ISBN: 0-8493-2479-3.
84
85 Parameters:
86
87 Input, int ORDER, the order.
88 ORDER must be between 2 and 20.
89
90 Output, double XTAB[ORDER], the abscissas.
91
92 Output, double WEIGHT[ORDER], the weights.
93 */
94 {
95 xtab.resize(order);
96 weight.resize(order);
97
98 if ( order == 2 )
99 {
100 xtab[0] = - 1.0E+00;
101 xtab[1] = 1.0E+00;
102
103 weight[0] = 1.0E+00;
104 weight[1] = 1.0E+00;
105 }
106 else if ( order == 3 )
107 {
108 xtab[0] = - 1.0E+00;
109 xtab[1] = 0.0E+00;
110 xtab[2] = 1.0E+00;
111
112 weight[0] = 1.0 / 3.0E+00;
113 weight[1] = 4.0 / 3.0E+00;
114 weight[2] = 1.0 / 3.0E+00;
115 }
116 else if ( order == 4 )
117 {
118 xtab[0] = - 1.0E+00;
119 xtab[1] = - 0.447213595499957939281834733746E+00;
120 xtab[2] = 0.447213595499957939281834733746E+00;
121 xtab[3] = 1.0E+00;
122
123 weight[0] = 1.0E+00 / 6.0E+00;
124 weight[1] = 5.0E+00 / 6.0E+00;
125 weight[2] = 5.0E+00 / 6.0E+00;
126 weight[3] = 1.0E+00 / 6.0E+00;
127 }
128 else if ( order == 5 )
129 {
130 xtab[0] = - 1.0E+00;
131 xtab[1] = - 0.654653670707977143798292456247E+00;
132 xtab[2] = 0.0E+00;
133 xtab[3] = 0.654653670707977143798292456247E+00;
134 xtab[4] = 1.0E+00;
135
136 weight[0] = 9.0E+00 / 90.0E+00;
137 weight[1] = 49.0E+00 / 90.0E+00;
138 weight[2] = 64.0E+00 / 90.0E+00;
139 weight[3] = 49.0E+00 / 90.0E+00;
140 weight[4] = 9.0E+00 / 90.0E+00;
141 }
142 else if ( order == 6 )
143 {
144 xtab[0] = - 1.0E+00;
145 xtab[1] = - 0.765055323929464692851002973959E+00;
146 xtab[2] = - 0.285231516480645096314150994041E+00;
147 xtab[3] = 0.285231516480645096314150994041E+00;
148 xtab[4] = 0.765055323929464692851002973959E+00;
149 xtab[5] = 1.0E+00;
150
151 weight[0] = 0.066666666666666666666666666667E+00;
152 weight[1] = 0.378474956297846980316612808212E+00;
153 weight[2] = 0.554858377035486353016720525121E+00;
154 weight[3] = 0.554858377035486353016720525121E+00;
155 weight[4] = 0.378474956297846980316612808212E+00;
156 weight[5] = 0.066666666666666666666666666667E+00;
157 }
158 else if ( order == 7 )
159 {
160 xtab[0] = - 1.0E+00;
161 xtab[1] = - 0.830223896278566929872032213967E+00;
162 xtab[2] = - 0.468848793470714213803771881909E+00;
163 xtab[3] = 0.0E+00;
164 xtab[4] = 0.468848793470714213803771881909E+00;
165 xtab[5] = 0.830223896278566929872032213967E+00;
166 xtab[6] = 1.0E+00;
167
168 weight[0] = 0.476190476190476190476190476190E-01;
169 weight[1] = 0.276826047361565948010700406290E+00;
170 weight[2] = 0.431745381209862623417871022281E+00;
171 weight[3] = 0.487619047619047619047619047619E+00;
172 weight[4] = 0.431745381209862623417871022281E+00;
173 weight[5] = 0.276826047361565948010700406290E+00;
174 weight[6] = 0.476190476190476190476190476190E-01;
175 }
176 else if ( order == 8 )
177 {
178 xtab[0] = - 1.0E+00;
179 xtab[1] = - 0.871740148509606615337445761221E+00;
180 xtab[2] = - 0.591700181433142302144510731398E+00;
181 xtab[3] = - 0.209299217902478868768657260345E+00;
182 xtab[4] = 0.209299217902478868768657260345E+00;
183 xtab[5] = 0.591700181433142302144510731398E+00;
184 xtab[6] = 0.871740148509606615337445761221E+00;
185 xtab[7] = 1.0E+00;
186
187 weight[0] = 0.357142857142857142857142857143E-01;
188 weight[1] = 0.210704227143506039382991065776E+00;
189 weight[2] = 0.341122692483504364764240677108E+00;
190 weight[3] = 0.412458794658703881567052971402E+00;
191 weight[4] = 0.412458794658703881567052971402E+00;
192 weight[5] = 0.341122692483504364764240677108E+00;
193 weight[6] = 0.210704227143506039382991065776E+00;
194 weight[7] = 0.357142857142857142857142857143E-01;
195 }
196 else if ( order == 9 )
197 {
198 xtab[0] = - 1.0E+00;
199 xtab[1] = - 0.899757995411460157312345244418E+00;
200 xtab[2] = - 0.677186279510737753445885427091E+00;
201 xtab[3] = - 0.363117463826178158710752068709E+00;
202 xtab[4] = 0.0E+00;
203 xtab[5] = 0.363117463826178158710752068709E+00;
204 xtab[6] = 0.677186279510737753445885427091E+00;
205 xtab[7] = 0.899757995411460157312345244418E+00;
206 xtab[8] = 1.0E+00;
207
208 weight[0] = 0.277777777777777777777777777778E-01;
209 weight[1] = 0.165495361560805525046339720029E+00;
210 weight[2] = 0.274538712500161735280705618579E+00;
211 weight[3] = 0.346428510973046345115131532140E+00;
212 weight[4] = 0.371519274376417233560090702948E+00;
213 weight[5] = 0.346428510973046345115131532140E+00;
214 weight[6] = 0.274538712500161735280705618579E+00;
215 weight[7] = 0.165495361560805525046339720029E+00;
216 weight[8] = 0.277777777777777777777777777778E-01;
217 }
218 else if ( order == 10 )
219 {
220 xtab[0] = - 1.0E+00;
221 xtab[1] = - 0.919533908166458813828932660822E+00;
222 xtab[2] = - 0.738773865105505075003106174860E+00;
223 xtab[3] = - 0.477924949810444495661175092731E+00;
224 xtab[4] = - 0.165278957666387024626219765958E+00;
225 xtab[5] = 0.165278957666387024626219765958E+00;
226 xtab[6] = 0.477924949810444495661175092731E+00;
227 xtab[7] = 0.738773865105505075003106174860E+00;
228 xtab[8] = 0.919533908166458813828932660822E+00;
229 xtab[9] = 1.0E+00;
230
231 weight[0] = 0.222222222222222222222222222222E-01;
232 weight[1] = 0.133305990851070111126227170755E+00;
233 weight[2] = 0.224889342063126452119457821731E+00;
234 weight[3] = 0.292042683679683757875582257374E+00;
235 weight[4] = 0.327539761183897456656510527917E+00;
236 weight[5] = 0.327539761183897456656510527917E+00;
237 weight[6] = 0.292042683679683757875582257374E+00;
238 weight[7] = 0.224889342063126452119457821731E+00;
239 weight[8] = 0.133305990851070111126227170755E+00;
240 weight[9] = 0.222222222222222222222222222222E-01;
241 }
242 else if ( order == 11 )
243 {
244 xtab[0] = - 1.0E+00;
245 xtab[1] = - 0.934001430408059134332274136099E+00;
246 xtab[2] = - 0.784483473663144418622417816108E+00;
247 xtab[3] = - 0.565235326996205006470963969478E+00;
248 xtab[4] = - 0.295758135586939391431911515559E+00;
249 xtab[5] = 0.0E+00;
250 xtab[6] = 0.295758135586939391431911515559E+00;
251 xtab[7] = 0.565235326996205006470963969478E+00;
252 xtab[8] = 0.784483473663144418622417816108E+00;
253 xtab[9] = 0.934001430408059134332274136099E+00;
254 xtab[10] = 1.0E+00;
255
256 weight[0] = 0.181818181818181818181818181818E-01;
257 weight[1] = 0.109612273266994864461403449580E+00;
258 weight[2] = 0.187169881780305204108141521899E+00;
259 weight[3] = 0.248048104264028314040084866422E+00;
260 weight[4] = 0.286879124779008088679222403332E+00;
261 weight[5] = 0.300217595455690693785931881170E+00;
262 weight[6] = 0.286879124779008088679222403332E+00;
263 weight[7] = 0.248048104264028314040084866422E+00;
264 weight[8] = 0.187169881780305204108141521899E+00;
265 weight[9] = 0.109612273266994864461403449580E+00;
266 weight[10] = 0.181818181818181818181818181818E-01;
267 }
268 else if ( order == 12 )
269 {
270 xtab[0] = - 1.0E+00;
271 xtab[1] = - 0.944899272222882223407580138303E+00;
272 xtab[2] = - 0.819279321644006678348641581717E+00;
273 xtab[3] = - 0.632876153031869677662404854444E+00;
274 xtab[4] = - 0.399530940965348932264349791567E+00;
275 xtab[5] = - 0.136552932854927554864061855740E+00;
276 xtab[6] = 0.136552932854927554864061855740E+00;
277 xtab[7] = 0.399530940965348932264349791567E+00;
278 xtab[8] = 0.632876153031869677662404854444E+00;
279 xtab[9] = 0.819279321644006678348641581717E+00;
280 xtab[10] = 0.944899272222882223407580138303E+00;
281 xtab[11] = 1.0E+00;
282
283 weight[0] = 0.151515151515151515151515151515E-01;
284 weight[1] = 0.916845174131961306683425941341E-01;
285 weight[2] = 0.157974705564370115164671062700E+00;
286 weight[3] = 0.212508417761021145358302077367E+00;
287 weight[4] = 0.251275603199201280293244412148E+00;
288 weight[5] = 0.271405240910696177000288338500E+00;
289 weight[6] = 0.271405240910696177000288338500E+00;
290 weight[7] = 0.251275603199201280293244412148E+00;
291 weight[8] = 0.212508417761021145358302077367E+00;
292 weight[9] = 0.157974705564370115164671062700E+00;
293 weight[10] = 0.916845174131961306683425941341E-01;
294 weight[11] = 0.151515151515151515151515151515E-01;
295 }
296 else if ( order == 13 )
297 {
298 xtab[0] = - 1.0E+00;
299 xtab[1] = - 0.953309846642163911896905464755E+00;
300 xtab[2] = - 0.846347564651872316865925607099E+00;
301 xtab[3] = - 0.686188469081757426072759039566E+00;
302 xtab[4] = - 0.482909821091336201746937233637E+00;
303 xtab[5] = - 0.249286930106239992568673700374E+00;
304 xtab[6] = 0.0E+00;
305 xtab[7] = 0.249286930106239992568673700374E+00;
306 xtab[8] = 0.482909821091336201746937233637E+00;
307 xtab[9] = 0.686188469081757426072759039566E+00;
308 xtab[10] = 0.846347564651872316865925607099E+00;
309 xtab[11] = 0.953309846642163911896905464755E+00;
310 xtab[12] = 1.0E+00;
311
312 weight[0] = 0.128205128205128205128205128205E-01;
313 weight[1] = 0.778016867468189277935889883331E-01;
314 weight[2] = 0.134981926689608349119914762589E+00;
315 weight[3] = 0.183646865203550092007494258747E+00;
316 weight[4] = 0.220767793566110086085534008379E+00;
317 weight[5] = 0.244015790306676356458578148360E+00;
318 weight[6] = 0.251930849333446736044138641541E+00;
319 weight[7] = 0.244015790306676356458578148360E+00;
320 weight[8] = 0.220767793566110086085534008379E+00;
321 weight[9] = 0.183646865203550092007494258747E+00;
322 weight[10] = 0.134981926689608349119914762589E+00;
323 weight[11] = 0.778016867468189277935889883331E-01;
324 weight[12] = 0.128205128205128205128205128205E-01;
325 }
326 else if ( order == 14 )
327 {
328 xtab[0] = - 1.0E+00;
329 xtab[1] = - 0.959935045267260901355100162015E+00;
330 xtab[2] = - 0.867801053830347251000220202908E+00;
331 xtab[3] = - 0.728868599091326140584672400521E+00;
332 xtab[4] = - 0.550639402928647055316622705859E+00;
333 xtab[5] = - 0.342724013342712845043903403642E+00;
334 xtab[6] = - 0.116331868883703867658776709736E+00;
335 xtab[7] = 0.116331868883703867658776709736E+00;
336 xtab[8] = 0.342724013342712845043903403642E+00;
337 xtab[9] = 0.550639402928647055316622705859E+00;
338 xtab[10] = 0.728868599091326140584672400521E+00;
339 xtab[11] = 0.867801053830347251000220202908E+00;
340 xtab[12] = 0.959935045267260901355100162015E+00;
341 xtab[13] = 1.0E+00;
342
343 weight[0] = 0.109890109890109890109890109890E-01;
344 weight[1] = 0.668372844976812846340706607461E-01;
345 weight[2] = 0.116586655898711651540996670655E+00;
346 weight[3] = 0.160021851762952142412820997988E+00;
347 weight[4] = 0.194826149373416118640331778376E+00;
348 weight[5] = 0.219126253009770754871162523954E+00;
349 weight[6] = 0.231612794468457058889628357293E+00;
350 weight[7] = 0.231612794468457058889628357293E+00;
351 weight[8] = 0.219126253009770754871162523954E+00;
352 weight[9] = 0.194826149373416118640331778376E+00;
353 weight[10] = 0.160021851762952142412820997988E+00;
354 weight[11] = 0.116586655898711651540996670655E+00;
355 weight[12] = 0.668372844976812846340706607461E-01;
356 weight[13] = 0.109890109890109890109890109890E-01;
357 }
358 else if ( order == 15 )
359 {
360 xtab[0] = - 1.0E+00;
361 xtab[1] = - 0.965245926503838572795851392070E+00;
362 xtab[2] = - 0.885082044222976298825401631482E+00;
363 xtab[3] = - 0.763519689951815200704118475976E+00;
364 xtab[4] = - 0.606253205469845711123529938637E+00;
365 xtab[5] = - 0.420638054713672480921896938739E+00;
366 xtab[6] = - 0.215353955363794238225679446273E+00;
367 xtab[7] = 0.0E+00;
368 xtab[8] = 0.215353955363794238225679446273E+00;
369 xtab[9] = 0.420638054713672480921896938739E+00;
370 xtab[10] = 0.606253205469845711123529938637E+00;
371 xtab[11] = 0.763519689951815200704118475976E+00;
372 xtab[12] = 0.885082044222976298825401631482E+00;
373 xtab[13] = 0.965245926503838572795851392070E+00;
374 xtab[14] = 1.0E+00;
375
376 weight[0] = 0.952380952380952380952380952381E-02;
377 weight[1] = 0.580298930286012490968805840253E-01;
378 weight[2] = 0.101660070325718067603666170789E+00;
379 weight[3] = 0.140511699802428109460446805644E+00;
380 weight[4] = 0.172789647253600949052077099408E+00;
381 weight[5] = 0.196987235964613356092500346507E+00;
382 weight[6] = 0.211973585926820920127430076977E+00;
383 weight[7] = 0.217048116348815649514950214251E+00;
384 weight[8] = 0.211973585926820920127430076977E+00;
385 weight[9] = 0.196987235964613356092500346507E+00;
386 weight[10] = 0.172789647253600949052077099408E+00;
387 weight[11] = 0.140511699802428109460446805644E+00;
388 weight[12] = 0.101660070325718067603666170789E+00;
389 weight[13] = 0.580298930286012490968805840253E-01;
390 weight[14] = 0.952380952380952380952380952381E-02;
391 }
392 else if ( order == 16 )
393 {
394 xtab[0] = - 1.0E+00;
395 xtab[1] = - 0.969568046270217932952242738367E+00;
396 xtab[2] = - 0.899200533093472092994628261520E+00;
397 xtab[3] = - 0.792008291861815063931088270963E+00;
398 xtab[4] = - 0.652388702882493089467883219641E+00;
399 xtab[5] = - 0.486059421887137611781890785847E+00;
400 xtab[6] = - 0.299830468900763208098353454722E+00;
401 xtab[7] = - 0.101326273521949447843033005046E+00;
402 xtab[8] = 0.101326273521949447843033005046E+00;
403 xtab[9] = 0.299830468900763208098353454722E+00;
404 xtab[10] = 0.486059421887137611781890785847E+00;
405 xtab[11] = 0.652388702882493089467883219641E+00;
406 xtab[12] = 0.792008291861815063931088270963E+00;
407 xtab[13] = 0.899200533093472092994628261520E+00;
408 xtab[14] = 0.969568046270217932952242738367E+00;
409 xtab[15] = 1.0E+00;
410
411 weight[0] = 0.833333333333333333333333333333E-02;
412 weight[1] = 0.508503610059199054032449195655E-01;
413 weight[2] = 0.893936973259308009910520801661E-01;
414 weight[3] = 0.124255382132514098349536332657E+00;
415 weight[4] = 0.154026980807164280815644940485E+00;
416 weight[5] = 0.177491913391704125301075669528E+00;
417 weight[6] = 0.193690023825203584316913598854E+00;
418 weight[7] = 0.201958308178229871489199125411E+00;
419 weight[8] = 0.201958308178229871489199125411E+00;
420 weight[9] = 0.193690023825203584316913598854E+00;
421 weight[10] = 0.177491913391704125301075669528E+00;
422 weight[11] = 0.154026980807164280815644940485E+00;
423 weight[12] = 0.124255382132514098349536332657E+00;
424 weight[13] = 0.893936973259308009910520801661E-01;
425 weight[14] = 0.508503610059199054032449195655E-01;
426 weight[15] = 0.833333333333333333333333333333E-02;
427 }
428 else if ( order == 17 )
429 {
430 xtab[0] = - 1.0E+00;
431 xtab[1] = - 0.973132176631418314156979501874E+00;
432 xtab[2] = - 0.910879995915573595623802506398E+00;
433 xtab[3] = - 0.815696251221770307106750553238E+00;
434 xtab[4] = - 0.691028980627684705394919357372E+00;
435 xtab[5] = - 0.541385399330101539123733407504E+00;
436 xtab[6] = - 0.372174433565477041907234680735E+00;
437 xtab[7] = - 0.189511973518317388304263014753E+00;
438 xtab[8] = 0.0E+00;
439 xtab[9] = 0.189511973518317388304263014753E+00;
440 xtab[10] = 0.372174433565477041907234680735E+00;
441 xtab[11] = 0.541385399330101539123733407504E+00;
442 xtab[12] = 0.691028980627684705394919357372E+00;
443 xtab[13] = 0.815696251221770307106750553238E+00;
444 xtab[14] = 0.910879995915573595623802506398E+00;
445 xtab[15] = 0.973132176631418314156979501874E+00;
446 xtab[16] = 1.0E+00;
447
448 weight[0] = 0.735294117647058823529411764706E-02;
449 weight[1] = 0.449219405432542096474009546232E-01;
450 weight[2] = 0.791982705036871191902644299528E-01;
451 weight[3] = 0.110592909007028161375772705220E+00;
452 weight[4] = 0.137987746201926559056201574954E+00;
453 weight[5] = 0.160394661997621539516328365865E+00;
454 weight[6] = 0.177004253515657870436945745363E+00;
455 weight[7] = 0.187216339677619235892088482861E+00;
456 weight[8] = 0.190661874753469433299407247028E+00;
457 weight[9] = 0.187216339677619235892088482861E+00;
458 weight[10] = 0.177004253515657870436945745363E+00;
459 weight[11] = 0.160394661997621539516328365865E+00;
460 weight[12] = 0.137987746201926559056201574954E+00;
461 weight[13] = 0.110592909007028161375772705220E+00;
462 weight[14] = 0.791982705036871191902644299528E-01;
463 weight[15] = 0.449219405432542096474009546232E-01;
464 weight[16] = 0.735294117647058823529411764706E-02;
465 }
466 else if ( order == 18 )
467 {
468 xtab[0] = - 1.0E+00;
469 xtab[1] = - 0.976105557412198542864518924342E+00;
470 xtab[2] = - 0.920649185347533873837854625431E+00;
471 xtab[3] = - 0.835593535218090213713646362328E+00;
472 xtab[4] = - 0.723679329283242681306210365302E+00;
473 xtab[5] = - 0.588504834318661761173535893194E+00;
474 xtab[6] = - 0.434415036912123975342287136741E+00;
475 xtab[7] = - 0.266362652878280984167665332026E+00;
476 xtab[8] = - 0.897490934846521110226450100886E-01;
477 xtab[9] = 0.897490934846521110226450100886E-01;
478 xtab[10] = 0.266362652878280984167665332026E+00;
479 xtab[11] = 0.434415036912123975342287136741E+00;
480 xtab[12] = 0.588504834318661761173535893194E+00;
481 xtab[13] = 0.723679329283242681306210365302E+00;
482 xtab[14] = 0.835593535218090213713646362328E+00;
483 xtab[15] = 0.920649185347533873837854625431E+00;
484 xtab[16] = 0.976105557412198542864518924342E+00;
485 xtab[17] = 1.0E+00;
486
487 weight[0] = 0.653594771241830065359477124183E-02;
488 weight[1] = 0.399706288109140661375991764101E-01;
489 weight[2] = 0.706371668856336649992229601678E-01;
490 weight[3] = 0.990162717175028023944236053187E-01;
491 weight[4] = 0.124210533132967100263396358897E+00;
492 weight[5] = 0.145411961573802267983003210494E+00;
493 weight[6] = 0.161939517237602489264326706700E+00;
494 weight[7] = 0.173262109489456226010614403827E+00;
495 weight[8] = 0.179015863439703082293818806944E+00;
496 weight[9] = 0.179015863439703082293818806944E+00;
497 weight[10] = 0.173262109489456226010614403827E+00;
498 weight[11] = 0.161939517237602489264326706700E+00;
499 weight[12] = 0.145411961573802267983003210494E+00;
500 weight[13] = 0.124210533132967100263396358897E+00;
501 weight[14] = 0.990162717175028023944236053187E-01;
502 weight[15] = 0.706371668856336649992229601678E-01;
503 weight[16] = 0.399706288109140661375991764101E-01;
504 weight[17] = 0.653594771241830065359477124183E-02;
505 }
506 else if ( order == 19 )
507 {
508 xtab[0] = - 1.0E+00;
509 xtab[1] = - 0.978611766222080095152634063110E+00;
510 xtab[2] = - 0.928901528152586243717940258797E+00;
511 xtab[3] = - 0.852460577796646093085955970041E+00;
512 xtab[4] = - 0.751494202552613014163637489634E+00;
513 xtab[5] = - 0.628908137265220497766832306229E+00;
514 xtab[6] = - 0.488229285680713502777909637625E+00;
515 xtab[7] = - 0.333504847824498610298500103845E+00;
516 xtab[8] = - 0.169186023409281571375154153445E+00;
517 xtab[9] = 0.0E+00;
518 xtab[10] = 0.169186023409281571375154153445E+00;
519 xtab[11] = 0.333504847824498610298500103845E+00;
520 xtab[12] = 0.488229285680713502777909637625E+00;
521 xtab[13] = 0.628908137265220497766832306229E+00;
522 xtab[14] = 0.751494202552613014163637489634E+00;
523 xtab[15] = 0.852460577796646093085955970041E+00;
524 xtab[16] = 0.928901528152586243717940258797E+00;
525 xtab[17] = 0.978611766222080095152634063110E+00;
526 xtab[18] = 1.0E+00;
527
528 weight[0] = 0.584795321637426900584795321637E-02;
529 weight[1] = 0.357933651861764771154255690351E-01;
530 weight[2] = 0.633818917626297368516956904183E-01;
531 weight[3] = 0.891317570992070844480087905562E-01;
532 weight[4] = 0.112315341477305044070910015464E+00;
533 weight[5] = 0.132267280448750776926046733910E+00;
534 weight[6] = 0.148413942595938885009680643668E+00;
535 weight[7] = 0.160290924044061241979910968184E+00;
536 weight[8] = 0.167556584527142867270137277740E+00;
537 weight[9] = 0.170001919284827234644672715617E+00;
538 weight[10] = 0.167556584527142867270137277740E+00;
539 weight[11] = 0.160290924044061241979910968184E+00;
540 weight[12] = 0.148413942595938885009680643668E+00;
541 weight[13] = 0.132267280448750776926046733910E+00;
542 weight[14] = 0.112315341477305044070910015464E+00;
543 weight[15] = 0.891317570992070844480087905562E-01;
544 weight[16] = 0.633818917626297368516956904183E-01;
545 weight[17] = 0.357933651861764771154255690351E-01;
546 weight[18] = 0.584795321637426900584795321637E-02;
547 }
548 else if ( order == 20 )
549 {
550 xtab[0] = - 1.0E+00;
551 xtab[1] = - 0.980743704893914171925446438584E+00;
552 xtab[2] = - 0.935934498812665435716181584931E+00;
553 xtab[3] = - 0.866877978089950141309847214616E+00;
554 xtab[4] = - 0.775368260952055870414317527595E+00;
555 xtab[5] = - 0.663776402290311289846403322971E+00;
556 xtab[6] = - 0.534992864031886261648135961829E+00;
557 xtab[7] = - 0.392353183713909299386474703816E+00;
558 xtab[8] = - 0.239551705922986495182401356927E+00;
559 xtab[9] = - 0.805459372388218379759445181596E-01;
560 xtab[10] = 0.805459372388218379759445181596E-01;
561 xtab[11] = 0.239551705922986495182401356927E+00;
562 xtab[12] = 0.392353183713909299386474703816E+00;
563 xtab[13] = 0.534992864031886261648135961829E+00;
564 xtab[14] = 0.663776402290311289846403322971E+00;
565 xtab[15] = 0.775368260952055870414317527595E+00;
566 xtab[16] = 0.866877978089950141309847214616E+00;
567 xtab[17] = 0.935934498812665435716181584931E+00;
568 xtab[18] = 0.980743704893914171925446438584E+00;
569 xtab[19] = 1.0E+00;
570
571 weight[0] = 0.526315789473684210526315789474E-02;
572 weight[1] = 0.322371231884889414916050281173E-01;
573 weight[2] = 0.571818021275668260047536271732E-01;
574 weight[3] = 0.806317639961196031447768461137E-01;
575 weight[4] = 0.101991499699450815683781205733E+00;
576 weight[5] = 0.120709227628674725099429705002E+00;
577 weight[6] = 0.136300482358724184489780792989E+00;
578 weight[7] = 0.148361554070916825814713013734E+00;
579 weight[8] = 0.156580102647475487158169896794E+00;
580 weight[9] = 0.160743286387845749007726726449E+00;
581 weight[10] = 0.160743286387845749007726726449E+00;
582 weight[11] = 0.156580102647475487158169896794E+00;
583 weight[12] = 0.148361554070916825814713013734E+00;
584 weight[13] = 0.136300482358724184489780792989E+00;
585 weight[14] = 0.120709227628674725099429705002E+00;
586 weight[15] = 0.101991499699450815683781205733E+00;
587 weight[16] = 0.806317639961196031447768461137E-01;
588 weight[17] = 0.571818021275668260047536271732E-01;
589 weight[18] = 0.322371231884889414916050281173E-01;
590 weight[19] = 0.526315789473684210526315789474E-02;
591 }
592 else
593 {
594 ERROR_INFO();
595 throw std::domain_error("Legal values for lobatto_set are between 2 and 20.\n");
596 }
597 }
598
599
lobatto_compute(int n,std::vector<double> & x,std::vector<double> & w)600 void lobatto_compute (int n, std::vector<double> & x, std::vector<double> & w)
601
602 /******************************************************************************/
603 /*
604 Purpose:
605
606 LOBATTO_COMPUTE computes a Lobatto quadrature rule.
607
608 Discussion:
609
610 The integral:
611
612 Integral ( -1 <= X <= 1 ) F(X) dX
613
614 The quadrature rule:
615
616 Sum ( 1 <= I <= N ) WEIGHT(I) * F ( XTAB(I) )
617
618 The quadrature rule will integrate exactly all polynomials up to
619 X**(2*N-3).
620
621 The Lobatto rule is distinguished by the fact that both endpoints
622 (-1 and 1) are always abscissas.
623
624 Licensing:
625
626 This code is distributed under the GNU LGPL license.
627
628 Modified:
629
630 04 February 2007
631
632 Author:
633
634 Original MATLAB version by Greg von Winckel.
635 C version by John Burkardt.
636
637 Reference:
638
639 Milton Abramowitz, Irene Stegun,
640 Handbook of Mathematical Functions,
641 National Bureau of Standards, 1964,
642 ISBN: 0-486-61272-4,
643 LC: QA47.A34.
644
645 Claudio Canuto, Yousuff Hussaini, Alfio Quarteroni, Thomas Zang,
646 Spectral Methods in Fluid Dynamics,
647 Springer, 1993,
648 ISNB13: 978-3540522058,
649 LC: QA377.S676.
650
651 Arthur Stroud, Don Secrest,
652 Gaussian Quadrature Formulas,
653 Prentice Hall, 1966,
654 LC: QA299.4G3S7.
655
656 Daniel Zwillinger, editor,
657 CRC Standard Mathematical Tables and Formulae,
658 30th Edition,
659 CRC Press, 1996,
660 ISBN: 0-8493-2479-3.
661
662 Parameters:
663
664 Input, int N, the order.
665 N must be at least 2.
666
667 Output, double X[N], the abscissas.
668
669 Output, double W[N], the weights.
670 */
671 {
672 int i;
673 int j;
674 double test, error;
675 double tolerance;
676
677 if ( n < 2 )
678 {
679 ERROR_INFO();
680 std::ostringstream oss;
681 oss << "Lobatto called with n="<<n<<", but n>=2 is required.\n";
682 throw std::runtime_error(oss.str());
683 }
684
685 if(n<20) {
686 // Use tabled weights and nodes
687 lobatto_set(n,x,w);
688 return;
689 }
690
691 // Resize tables to correct length
692 x.resize(n);
693 w.resize(n);
694
695 tolerance = 100.0 * DBL_EPSILON;
696 /*
697 Initial estimate for the abscissas is the Chebyshev-Gauss-Lobatto nodes.
698 */
699 for ( i = 0; i < n; i++ )
700 {
701 x[i] = cos ( M_PI * ( double ) ( i ) / ( double ) ( n - 1 ) );
702 }
703
704 double xold[n];
705 double p[n*n];
706
707 do
708 {
709 for ( i = 0; i < n; i++ )
710 {
711 xold[i] = x[i];
712 }
713 for ( i = 0; i < n; i++ )
714 {
715 p[i+0*n] = 1.0;
716 }
717 for ( i = 0; i < n; i++ )
718 {
719 p[i+1*n] = x[i];
720 }
721
722 for ( j = 2; j <= n-1; j++ )
723 {
724 for ( i = 0; i < n; i++)
725 {
726 p[i+j*n] = ( ( double ) ( 2 * j - 1 ) * x[i] * p[i+(j-1)*n]
727 + ( double ) ( - j + 1 ) * p[i+(j-2)*n] )
728 / ( double ) ( j );
729 }
730 }
731
732 for ( i = 0; i < n; i++ )
733 {
734 x[i] = xold[i] - ( x[i] * p[i+(n-1)*n] - p[i+(n-2)*n] )
735 / ( ( double ) ( n ) * p[i+(n-1)*n] );
736 }
737
738 error = 0.0;
739 for ( i = 0; i < n; i++ )
740 {
741 test = fabs(x[i] - xold[i]);
742 if(test>error)
743 error=test;
744 }
745
746 } while ( tolerance < error );
747
748 // Reverse order of x
749 reverse(x);
750
751 for ( i = 0; i < n; i++ )
752 {
753 w[i] = 2.0 / ( ( double ) ( ( n - 1 ) * n ) * pow ( p[i+(n-1)*n], 2 ) );
754 }
755 }
756 /******************************************************************************/
757
758