1 /*
2 XLiFE++ is an extended library of finite elements written in C++
3     Copyright (C) 2014  Lunéville, Eric; Kielbasiewicz, Nicolas; Lafranche, Yvon; Nguyen, Manh-Ha; Chambeyron, Colin
4 
5     This program is free software: you can redistribute it and/or modify
6     it under the terms of the GNU General Public License as published by
7     the Free Software Foundation, either version 3 of the License, or
8     (at your option) any later version.
9     This program is distributed in the hope that it will be useful,
10     but WITHOUT ANY WARRANTY; without even the implied warranty of
11     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12     GNU General Public License for more details.
13     You should have received a copy of the GNU General Public License
14     along with this program.  If not, see <http://www.gnu.org/licenses/>.
15  */
16 
17 /*!
18   \file GaussFormulae.cpp
19   \author D. Martin
20   \since 10 jan 2003
21   \date 4 jul 2011
22 
23   \brief Implementation of Gauss formulae functions
24  */
25 
26 #include "GaussFormulae.hpp"
27 
28 namespace xlifepp
29 {
30 
31 /*!
32   Compute Gauss-Legendre n-point formula on [-1, 1]
33   - Points xj are roots of Legendre Polynomial pn of ordre n
34   - Weights are 2/( (1-xj^2) (P'_n(xj)^2) )
35 
36   Note : the (n+1)/2 first positive points in ascending order only are returned.
37  */
gaussLegendreRuleComputed(number_t n,std::vector<real_t> & points,std::vector<real_t> & weights)38 void gaussLegendreRuleComputed(number_t n, std::vector<real_t>& points, std::vector<real_t>& weights)
39 {
40   const int nhalf((n + 1) / 2);
41   real_t xj(0.), x1(0.), delta;
42 
43   for (int j = 0; j < nhalf; j++)
44   {
45     real_t pn, pn1, pn2, ppn;
46     // Starting point of Newton algorithm for j-th root
47     xj = std::cos(pi_ * (j + 0.75) / (n + 0.5));
48     do
49     {
50       // compute Legendre polynomials P_{n-2}, P_{n-1}, P_{n} by recursion :
51       // k*P_k(x)=(2*k-1)*x*P_{k-1}(x) - (k-1)*P_{k-2}(x) , k=1..n
52       pn = 1.; pn1 = 0.;
53       for (number_t k = 1, km1 = 0, kx2m1 = 1; k <= n; k++, km1++, kx2m1 += 2)
54       {
55         pn2 = pn1;
56         pn1 = pn;
57         pn = (kx2m1 * xj * pn1 - km1 * pn2) / k;
58       }
59 
60       // compute derivative for Newton algorithm
61       //
62       // compute (1-x^2)*P'_n(x) at xj :
63       // (1-x^2) P'_n(x)=n * P_{n-1}(x) - n * x * pn(x)
64       ppn = n * (pn1 - xj * pn);
65       // Newton iteration update
66       delta = (1 - xj * xj) * pn / ppn;
67       x1 = xj;
68       xj -= delta;
69     }
70     while (std::abs(delta) > theEpsilon);
71     // previous solution is good enough (or else change tolerance)
72     weights[nhalf - 1 - j] = 2.*(1. - x1 * x1) / (ppn * ppn);
73     points[nhalf - 1 - j] = x1;
74   }
75   // Very small negative (near zero) values seem to occur for some odd n
76   // In case we want only non-negative values ...
77   // if ( abs(points[0]) < theEpsilon ) points[0]=0.
78   // As a matter of fact, we might as well set points[0]=0 for all odd n
79   // if ( n%2 == 1 ) points[0]=0;
80 }
81 
82 //! Output of Gauss-Legendre rule as displayed in function GaussLegendreRule
gaussLegendreOutput(number_t nmax,std::ostream & out)83 void gaussLegendreOutput(number_t nmax, std::ostream& out)
84 {
85   std::vector<real_t> p((nmax + 1) / 2), w((nmax + 1) / 2);
86   // Known exactly for n=1, 2, 3 and points only for n=4, 5
87   for (number_t n = 2; n < nmax; n++)
88   {
89     out << "      case " << n << ":" << std::endl;
90     gaussLegendreRule(n, p, w);
91     for (number_t i = 0; i < (n + 1) / 2; i++)
92     {
93       out.setf(std::ios::scientific);
94       out << "         *pi++=" << std::setprecision(19) << p[i] << "L; *wi=" << w[i] << "L;" << std::endl;
95     }
96     out << "         break;" << std::endl;
97     gaussLegendreRuleComputed(n, p, w);
98     for (number_t i = 0; i < (n + 1) / 2; i++)
99     {
100       out.setf(std::ios::scientific);
101       out << "         *pi++=" << std::setprecision(19) << p[i] << "L; *wi=" << w[i] << "L;" << std::endl;
102     }
103     out << "         break;// computed" << std::endl;
104   }
105 }
106 
107 /*!
108   returns quadrature points and weights for Gauss-Legendre formula for a given
109   number of points n ; n-point Gauss-Legendre formula is exact for
110   polynomials of degree up to {2n-1} on [-1,1].
111   Points and weights are given below for n up to 16
112 
113   Note : the (n+1)/2 positive points in ascending order only are returned.
114  */
gaussLegendreRule(number_t n,std::vector<real_t> & points,std::vector<real_t> & weights)115 void gaussLegendreRule(number_t n, std::vector<real_t>& points, std::vector<real_t>& weights)
116 {
117   // a few Legendre polynomials with easy roots
118   // P_2=(1/2) (3x^2-1);         P'_2=  3x
119   // P_3=(x/2) (5x^2-3)          P'_3= (3/2) (5x^2-1)
120   // P_4=(1/8) (35x^4-30x^2+3)   P'_4=(20x/8)(7x^2-3)
121   // P_5=(x/8) (63x^4-70x^2+15)  P'_5=(15/8) (21x^4-14x^2+1)
122   // weight associated with root r is 2.*(1.-r*r)/(P'_n(r))^2;
123   std::vector<real_t>::iterator it_p(points.begin()), it_w(weights.begin());
124   switch (n)
125   {
126     case 1: // 1 point { 0 } weight { 2 } (mid-point rule)
127       *it_p++ = 0.0L;                                 *it_w++ = 2.0L;
128       break;
129     case 2: // 2 points : { 1/sqrt3 }, weights { 1 }
130       *it_p++ = std::sqrt(real_t(1.) / 3);            *it_w++ = 1.0L;
131       // *it_p++=5.77350269189625764509148780501e-01L; *it_w++=1.0L;
132       break;
133     case 3: // 3 points {  0 sqrt{3/5} }, weights { 8/9  5/9 }
134       *it_p++ = 0.0L;                                 *it_w++ = real_t(8.) / 9; //8.8888888888888888888e-01L;
135       *it_p++ = std::sqrt(real_t(3.) / 5);            *it_w++ = real_t(5.) / 9; //5.5555555555555555555e-01L;
136       // *it_p++=7.74596669241483377035853079956e-01L; *it_w++=5.55555555555555555555555555555e-01L;
137       break;
138     case 4: // 4 points { sqrt((15[-/+]sqrt(120))/35) }, weights { }
139       *it_p++ = 3.39981043584856264802665759103e-01L; *it_w++ = 6.52145154862546142626936050778e-01L;
140       *it_p++ = 8.61136311594052575223946488893e-01L; *it_w++ = 3.47854845137453857373063949222e-01L;
141       break;
142     case 5: // 5 points { 0 sqrt((35[+/-]sqrt(280))/63) }, weights { }
143       *it_p++ = 0.L;                                  *it_w++ = 5.68888888888888888888888888889e-01L;
144       *it_p++ = 5.38469310105683091036314420700e-01L; *it_w++ = 4.78628670499366468041291514836e-01L;
145       *it_p++ = 9.06179845938663992797626878299e-01L; *it_w++ = 2.36926885056189087514264040720e-01L;
146       break;
147     case 6:
148       *it_p++ = 2.38619186083196908630501721681e-01L; *it_w++ = 4.67913934572691047389870343990e-01L;
149       *it_p++ = 6.61209386466264513661399595020e-01L; *it_w++ = 3.60761573048138607569833513838e-01L;
150       *it_p++ = 9.32469514203152027812301554494e-01L; *it_w++ = 1.71324492379170345040296142173e-01L;
151       break;
152     case 7:
153       *it_p++ = 0.L;                                  *it_w++ = 4.17959183673469387755102040816e-01L;
154       *it_p++ = 4.05845151377397166906606412077e-01L; *it_w++ = 3.81830050505118944950369775489e-01L;
155       *it_p++ = 7.41531185599394439863864773281e-01L; *it_w++ = 2.79705391489276667901467771424e-01L;
156       *it_p++ = 9.49107912342758524526189684048e-01L; *it_w++ = 1.29484966168869693270611432679e-01L;
157       break;
158     case 8:
159       *it_p++ = 1.83434642495649804939476142360e-01L; *it_w++ = 3.62683783378361982965150449277e-01L;
160       *it_p++ = 5.25532409916328985817739049189e-01L; *it_w++ = 3.13706645877887287337962201987e-01L;
161       *it_p++ = 7.96666477413626739591553936476e-01L; *it_w++ = 2.22381034453374470544355994426e-01L;
162       *it_p++ = 9.60289856497536231683560868569e-01L; *it_w++ = 1.01228536290376259152531354310e-01L;
163       break;
164     case 9:
165       *it_p++ = 0.L;                                  *it_w++ = 3.30239355001259763164525069287e-01L;
166       *it_p++ = 3.24253423403808929038538014643e-01L; *it_w++ = 3.12347077040002840068630406584e-01L;
167       *it_p++ = 6.13371432700590397308702039341e-01L; *it_w++ = 2.60610696402935462318742869419e-01L;
168       *it_p++ = 8.36031107326635794299429788070e-01L; *it_w++ = 1.80648160694857404058472031243e-01L;
169       *it_p++ = 9.68160239507626089835576202904e-01L; *it_w++ = 8.12743883615744119718921581105e-02L;
170       break;
171     case 10:
172       *it_p++ = 1.48874338981631210884826001130e-01L; *it_w++ = 2.95524224714752870173892994651e-01L;
173       *it_p++ = 4.33395394129247190799265943166e-01L; *it_w++ = 2.69266719309996355091226921569e-01L;
174       *it_p++ = 6.79409568299024406234327365115e-01L; *it_w++ = 2.19086362515982043995534934228e-01L;
175       *it_p++ = 8.65063366688984510732096688423e-01L; *it_w++ = 1.49451349150580593145776339658e-01L;
176       *it_p++ = 9.73906528517171720077964012084e-01L; *it_w++ = 6.66713443086881375935688098933e-02L;
177       break;
178     case 11:
179       *it_p++ = 0.L;                                  *it_w++ = 2.72925086777900630714483528336e-01L;
180       *it_p++ = 2.69543155952344972331531985401e-01L; *it_w++ = 2.62804544510246662180688869891e-01L;
181       *it_p++ = 5.19096129206811815925725669459e-01L; *it_w++ = 2.33193764591990479918523704843e-01L;
182       *it_p++ = 7.30152005574049324093416252031e-01L; *it_w++ = 1.86290210927734251426097641432e-01L;
183       *it_p++ = 8.87062599768095299075157769304e-01L; *it_w++ = 1.25580369464904624634694299224e-01L;
184       *it_p++ = 9.78228658146056992803938001123e-01L; *it_w++ = 5.56685671161736664827537204425e-02L;
185       break;
186     case 12:
187       *it_p++ = 1.25233408511468915472441369464e-01L; *it_w++ = 2.49147045813402785000562436043e-01L;
188       *it_p++ = 3.67831498998180193752691536644e-01L; *it_w++ = 2.33492536538354808760849898925e-01L;
189       *it_p++ = 5.87317954286617447296702418941e-01L; *it_w++ = 2.03167426723065921749064455810e-01L;
190       *it_p++ = 7.69902674194304687036893833213e-01L; *it_w++ = 1.60078328543346226334652529543e-01L;
191       *it_p++ = 9.04117256370474856678465866119e-01L; *it_w++ = 1.06939325995318430960254718194e-01L;
192       *it_p++ = 9.81560634246719250690549090149e-01L; *it_w++ = 4.71753363865118271946159614850e-02L;
193       break;
194     case 13:
195       *it_p++ = 0.L;                                  *it_w++ = 2.32551553230873910194589515269e-01L;
196       *it_p++ = 2.30458315955134794065528121098e-01L; *it_w++ = 2.26283180262897238412090186040e-01L;
197       *it_p++ = 4.48492751036446852877912852128e-01L; *it_w++ = 2.07816047536888502312523219306e-01L;
198       *it_p++ = 6.42349339440340220643984606996e-01L; *it_w++ = 1.78145980761945738280046691996e-01L;
199       *it_p++ = 8.01578090733309912794206489583e-01L; *it_w++ = 1.38873510219787238463601776869e-01L;
200       *it_p++ = 9.17598399222977965206547836501e-01L; *it_w++ = 9.21214998377284479144217759538e-02L;
201       *it_p++ = 9.84183054718588149472829448807e-01L; *it_w++ = 4.04840047653158795200215922010e-02L;
202       break;
203     case 14:
204       *it_p++ = 1.08054948707343662066244650220e-01L; *it_w++ = 2.15263853463157790195876443316e-01L;
205       *it_p++ = 3.19112368927889760435671824168e-01L; *it_w++ = 2.05198463721295603965924065661e-01L;
206       *it_p++ = 5.15248636358154091965290718551e-01L; *it_w++ = 1.85538397477937813741716590125e-01L;
207       *it_p++ = 6.87292904811685470148019803019e-01L; *it_w++ = 1.57203167158193534569601938624e-01L;
208       *it_p++ = 8.27201315069764993189794742650e-01L; *it_w++ = 1.21518570687903184689414809072e-01L;
209       *it_p++ = 9.28434883663573517336391139378e-01L; *it_w++ = 8.01580871597602098056332770629e-02L;
210       *it_p++ = 9.86283808696812338841597266704e-01L; *it_w++ = 3.51194603317518630318328761382e-02L;
211       break;
212     case 15:
213       *it_p++ = 0.L;                                  *it_w++ = 2.02578241925561272880620199968e-01L;
214       *it_p++ = 2.01194093997434522300628303395e-01L; *it_w++ = 1.98431485327111576456118326444e-01L;
215       *it_p++ = 3.94151347077563369897207370981e-01L; *it_w++ = 1.86161000015562211026800561866e-01L;
216       *it_p++ = 5.70972172608538847537226737254e-01L; *it_w++ = 1.66269205816993933553200860481e-01L;
217       *it_p++ = 7.24417731360170047416186054614e-01L; *it_w++ = 1.39570677926154314447804794511e-01L;
218       *it_p++ = 8.48206583410427216200648320774e-01L; *it_w++ = 1.07159220467171935011869546686e-01L;
219       *it_p++ = 9.37273392400705904307758947710e-01L; *it_w++ = 7.03660474881081247092674164507e-02L;
220       *it_p++ = 9.87992518020485428489565718587e-01L; *it_w++ = 3.07532419961172683546283935772e-02L;
221       break;
222     case 16:
223       *it_p++ = 9.50125098376374401853193354250e-02L; *it_w++ = 1.89450610455068496285396723208e-01L;
224       *it_p++ = 2.81603550779258913230460501460e-01L; *it_w++ = 1.82603415044923588866763667969e-01L;
225       *it_p++ = 4.58016777657227386342419442984e-01L; *it_w++ = 1.69156519395002538189312079030e-01L;
226       *it_p++ = 6.17876244402643748446671764049e-01L; *it_w++ = 1.49595988816576732081501730547e-01L;
227       *it_p++ = 7.55404408355003033895101194847e-01L; *it_w++ = 1.24628971255533872052476282192e-01L;
228       *it_p++ = 8.65631202387831743880467897712e-01L; *it_w++ = 9.51585116824927848099251076022e-02L;
229       *it_p++ = 9.44575023073232576077988415535e-01L; *it_w++ = 6.22535239386478928628438369944e-02L;
230       *it_p++ = 9.89400934991649932596154173450e-01L; *it_w++ = 2.71524594117540948517805724560e-02L;
231       break;
232     default:
233       gaussLegendreRuleComputed(n, points, weights);
234       break;
235   }
236 }
237 
238 /*!
239   compute Gauss-Lobatto n-point formula on [-1, 1] end points included
240     - Points xj are roots of (1-x^2)*P'_{n-1}
241                   with P'_{n-1} derivative of Legendre Polynomial of ordre n-1
242     - Weights are 2./(n*(n-1)) for end points -1 and +1
243                   2./(n*(n-1)*{P_{n-1}(xj)}^2) for other xj
244 
245    Note : the (n-1)/2 first positive points in ascending order only are returned
246         that is excluding end point 1 and corresponding weight.
247  */
gaussLobattoRuleComputed(number_t n,std::vector<real_t> & points,std::vector<real_t> & weights)248 void gaussLobattoRuleComputed(number_t n, std::vector<real_t>& points, std::vector<real_t>& weights)
249 {
250   points[0] = 0.; weights[0] = 2.;
251   const int nhalf((n + 1) / 2), nm1(n - 1), nnm1(n * nm1);
252   real_t xj(0.), x1(0.), delta(0);
253   int nj( std::max(0, nhalf - 2));
254   for (int j = 1; j < nhalf; j++, nj--)
255   {
256     real_t pn, pn1, pn2;
257     // Starting point of Newton algorithm for j-th root
258     xj = std::cos(pi_ * (j + 0.5) / (n + 0.5));
259     do
260     {
261       // compute Legendre polynomials P_{n-3}, P_{n-2}, P_{n-1} by recursion :
262       // k*P_k(x)=(2*k-1)*x*P_{k-1}(x) - (k-1)*P_{k-2}(x) , k=1..n-1
263       pn = 1.; pn1 = 0.;
264       for ( number_t k = 1, km1 = 0, kx2m1 = 1; k < n; k++, km1++, kx2m1 += 2)
265       {
266         pn2 = pn1;
267         pn1 = pn;
268         pn = (kx2m1 * xj * pn1 - km1 * pn2) / k;
269       }
270       // compute derivatives for Newton algorithm
271       // (1-x^2)*P'_{n-1}(x)=(n-1) * P_{n-2}(x) - (n-1) * x * P_{n-1}(x)
272       real_t ppn = nm1 * (pn1 - xj * pn);
273       // (1-x^2)*P''_{n-1}(x)=2 * x * P'_{n-1}(x) - n * (n-1) * P_{n-1}(x)
274       real_t ppn2 = 2.*xj * ppn / (1. - xj * xj) - nnm1 * pn;
275       // Newton iteration update
276       delta = ppn / ppn2;
277       x1 = xj;
278       xj -= delta;
279     }
280     while (std::abs(delta) > theEpsilon);
281     // previous solution is good enough (or else change tolerance)
282     weights[nj] = 2. / (nnm1 * pn * pn);
283     points[nj] = x1;
284   }
285   // Very small negative (near zero) values seem to occur for some odd n=19, 29, ...
286   // In case we want only non-negative values ...
287   // if ( abs(points[0]) < theEpsilon ) points[0]=0.;
288   // As a matter of fact, we might as well set points[0]=0 for all odd n
289   // if ( n%2 == 1 ) points[0]=0;
290 }
291 
292 //! Output of Gauss-Lobatto rule as displayed in function GaussLobattoRule
gaussLobattoOutput(number_t nmax,std::ostream & out)293 void gaussLobattoOutput(number_t nmax, std::ostream& out)
294 {
295   std::vector<real_t> p((nmax + 1) / 2), w((nmax + 1) / 2);
296   // Known exactly for n=2, 3, 4, 5
297   for (number_t n = 6; n < nmax; n++)
298   {
299     out << "      case " << n << ":" << std::endl;
300     gaussLobattoRule(n, p, w);
301     for ( number_t i = 0; i < (n + 1) / 2 - 1; i++)
302     {
303       out.setf(std::ios::scientific);
304       out << "         *pi=" << std::setprecision(19) << p[i] << "L; *wi=" << w[i] << "L;" << std::endl;
305     }
306     out << "         break;" << std::endl;
307     gaussLobattoRuleComputed(n, p, w);
308     for ( number_t i = 0; i < (n + 1) / 2 - 1; i++)
309     {
310       out.setf(std::ios::scientific);
311       out << "         *pi=" << std::setprecision(19) << p[i] << "L; *wi=" << w[i] << "L;" << std::endl;
312     }
313     out << "         break;// computed" << std::endl;
314   }
315 }
316 
317 /*!
318   returns quadrature points & weights for Gauss-Lobatto rule for a given number of points.
319   n-point Gauss-Lobatto formula is exact for P_{2n-3}(0,1)
320 
321   Note : only points and weights corresponding to the (n+1)/2 points non-negative points,
322   with possible point 0 included and including end point 1, are returned
323 
324   The following programs shows table of values (with 20 significant digits) up to n=16
325   which have been computed by function GaussLobattoRuleComputed using type long double.
326  */
gaussLobattoRule(number_t n,std::vector<real_t> & points,std::vector<real_t> & weights)327 void gaussLobattoRule(number_t n, std::vector<real_t>& points, std::vector<real_t>& weights)
328 {
329   // a few Legendre polynomials with easy roots
330   // P_2=(1/2) (3x^2-1);         P'_2=  3x
331   // P_3=(x/2) (5x^2-3)          P'_3= (3/2) (5x^2-1)
332   // P_4=(1/8) (35x^4-30x^2+3)   P'_4=(20x/8)(7x^2-3)
333   // P_5=(x/8) (63x^4-70x^2+15)  P'_5=(15/8) (21x^4-14x^2+1)
334   // weight associated with root r of P'_n is 2./( n*(n-1)*(pn(r))^2 )
335   std::vector<real_t>::iterator it_p(points.begin()), it_w(weights.begin());
336   switch (n)
337   {
338     case 2: // 2 point { 1 } (trapezoidal rule)
339       break;
340     case 3: // 3 point { 0  1 }, weights { 4/3 1/3 } (Simpson's rule)
341       *it_p++ = 0.0L;                       *it_w++ = real_t(4.) / 3; // 1.3333333333333333333e_01L;
342       break;
343     case 4: // 4 points { 1/sqrt5 1 }, weights { 5/6 1/6 }
344       *it_p++ = std::sqrt(real_t(1.) / 5);  *it_w++ = real_t(5.) / 6;
345       // *it_p++=4.47213595499957939281834733746e-01L; *it_w++=8.3333333333333333333e-01L;
346       break;
347     case 5: // 5 points { 0 sqrt(3/7) 1 }, weights { 32/45 49/90 1/10 }
348       *it_p++ = 0.0L;                       *it_w++ = real_t(32.) / 45; // 7.1111111111111111111e-01L;
349       *it_p++ = std::sqrt(real_t(3.) / 7);  *it_w++ = real_t(49.) / 90;
350       // *it_p++=6.54653670707977143798292456246e-01L; *it_w++=5.4444444444444444444e-01L;
351       break;
352     case 6:
353       *it_p++ = 2.8523151648064509631e-01L; *it_w++ = 5.5485837703548635306e-01L;
354       *it_p++ = 7.6505532392946469283e-01L; *it_w++ = 3.7847495629784698028e-01L;
355       break;
356     case 7:
357       *it_p++ = 0.0L;                       *it_w++ = 4.8761904761904761906e-01L;
358       *it_p++ = 4.6884879347071421380e-01L; *it_w++ = 4.3174538120986262336e-01L;
359       *it_p++ = 8.3022389627856692983e-01L; *it_w++ = 2.7682604736156594803e-01L;
360       break;
361     case 8:
362       *it_p++ = 2.0929921790247886877e-01L; *it_w++ = 4.1245879465870388160e-01L;
363       *it_p++ = 5.9170018143314230216e-01L; *it_w++ = 3.4112269248350436484e-01L;
364       *it_p++ = 8.7174014850960661533e-01L; *it_w++ = 2.1070422714350603931e-01L;
365       break;
366     case 9:
367       *it_p++ = 0.0L;                       *it_w++ = 3.7151927437641723355e-01L;
368       *it_p++ = 3.6311746382617815873e-01L; *it_w++ = 3.4642851097304634518e-01L;
369       *it_p++ = 6.7718627951073775344e-01L; *it_w++ = 2.7453871250016173528e-01L;
370       *it_p++ = 8.9975799541146015732e-01L; *it_w++ = 1.6549536156080552509e-01L;
371       break;
372     case 10:
373       *it_p++ = 1.6527895766638702462e-01L; *it_w++ = 3.2753976118389745660e-01L;
374       *it_p++ = 4.7792494981044449567e-01L; *it_w++ = 2.9204268367968375789e-01L;
375       *it_p++ = 7.3877386510550507501e-01L; *it_w++ = 2.2488934206312645217e-01L;
376       *it_p++ = 9.1953390816645881383e-01L; *it_w++ = 1.3330599085107011112e-01L;
377       break;
378     case 11:
379       *it_p++ = 0.0L;                       *it_w++ = 3.0021759545569069377e-01L;
380       *it_p++ = 2.9575813558693939142e-01L; *it_w++ = 2.8687912477900808858e-01L;
381       *it_p++ = 5.6523532699620500645e-01L; *it_w++ = 2.4804810426402831406e-01L;
382       *it_p++ = 7.8448347366314441861e-01L; *it_w++ = 1.8716988178030520414e-01L;
383       *it_p++ = 9.3400143040805913433e-01L; *it_w++ = 1.0961227326699486447e-01L;
384       break;
385     case 12:
386       *it_p++ = 1.3655293285492755487e-01L; *it_w++ = 2.7140524091069617696e-01L;
387       *it_p++ = 3.9953094096534893227e-01L; *it_w++ = 2.5127560319920128027e-01L;
388       *it_p++ = 6.3287615303186067767e-01L; *it_w++ = 2.1250841776102114536e-01L;
389       *it_p++ = 8.1927932164400667833e-01L; *it_w++ = 1.5797470556437011515e-01L;
390       *it_p++ = 9.4489927222288222340e-01L; *it_w++ = 9.1684517413196130668e-02L;
391       break;
392     case 13:
393       *it_p++ = 0.0L;                       *it_w++ = 2.5193084933344673604e-01L;
394       *it_p++ = 2.4928693010623999258e-01L; *it_w++ = 2.4401579030667635648e-01L;
395       *it_p++ = 4.8290982109133620174e-01L; *it_w++ = 2.2076779356611008604e-01L;
396       *it_p++ = 6.8618846908175742608e-01L; *it_w++ = 1.8364686520355009194e-01L;
397       *it_p++ = 8.4634756465187231688e-01L; *it_w++ = 1.3498192668960834912e-01L;
398       *it_p++ = 9.5330984664216391190e-01L; *it_w++ = 7.7801686746818927736e-02L;
399       break;
400     case 14:
401       *it_p++ = 1.1633186888370386766e-01L; *it_w++ = 2.3161279446845705891e-01L;
402       *it_p++ = 3.4272401334271284504e-01L; *it_w++ = 2.1912625300977075490e-01L;
403       *it_p++ = 5.5063940292864705533e-01L; *it_w++ = 1.9482614937341611863e-01L;
404       *it_p++ = 7.2886859909132614059e-01L; *it_w++ = 1.6002185176295214241e-01L;
405       *it_p++ = 8.6780105383034725100e-01L; *it_w++ = 1.1658665589871165151e-01L;
406       *it_p++ = 9.5993504526726090137e-01L; *it_w++ = 6.6837284497681284639e-02L;
407       break;
408     case 15:
409       *it_p++ = 0.0L;                       *it_w++ = 2.1704811634881564952e-01L;
410       *it_p++ = 2.1535395536379423822e-01L; *it_w++ = 2.1197358592682092013e-01L;
411       *it_p++ = 4.2063805471367248092e-01L; *it_w++ = 1.9698723596461335611e-01L;
412       *it_p++ = 6.0625320546984571109e-01L; *it_w++ = 1.7278964725360094905e-01L;
413       *it_p++ = 7.6351968995181520070e-01L; *it_w++ = 1.4051169980242810947e-01L;
414       *it_p++ = 8.8508204422297629880e-01L; *it_w++ = 1.0166007032571806761e-01L;
415       *it_p++ = 9.6524592650383857283e-01L; *it_w++ = 5.8029893028601249073e-02L;
416       break;
417     case 16:
418       *it_p++ = 1.0132627352194944784e-01L; *it_w++ = 2.0195830817822987145e-01L;
419       *it_p++ = 2.9983046890076320810e-01L; *it_w++ = 1.9369002382520358435e-01L;
420       *it_p++ = 4.8605942188713761179e-01L; *it_w++ = 1.7749191339170412529e-01L;
421       *it_p++ = 6.5238870288249308947e-01L; *it_w++ = 1.5402698080716428083e-01L;
422       *it_p++ = 7.9200829186181506393e-01L; *it_w++ = 1.2425538213251409831e-01L;
423       *it_p++ = 8.9920053309347209301e-01L; *it_w++ = 8.9393697325930800996e-02L;
424       *it_p++ = 9.6956804627021793294e-01L; *it_w++ = 5.0850361005919905360e-02L;
425       break;
426     default:
427       gaussLobattoRuleComputed(n, points, weights);
428       it_p = points.begin() + (n - 1) / 2;
429       it_w = weights.begin() + (n - 1) / 2;
430       break;
431   }
432   // end point ( abscissa=1, weight=2./(n(n-1)) )
433   *it_p = 1.L;
434   *it_w = 2. / (n * (n - 1));
435 }
436 
437 /*!
438   returns quadrature nodes for Gauss-Lobatto rule for a given number of points.
439 
440   Note : the (n+1)/2 points non-negative points in ascending order are returned only
441   with possible point 0 included and including end point 1
442  */
gaussLobattoPoints(number_t n,std::vector<real_t> & points)443 void gaussLobattoPoints(number_t n, std::vector<real_t>& points)
444 {
445   std::vector<real_t> weights((n + 1) / 2);
446   gaussLobattoRule(n, points, weights);
447 }
448 
449 /*!
450   Compute Gauss-Jacobi n-point formula on [-1, 1]
451   - Points xj are roots of Jacobi Polynomial Pn(a,b) of ordre n (namely orthogonal polynomials of weight (1-x)^a*(1+x)^b)
452   - Weights are - (2n+a+b+2)*Gamma(n+a+1)*Gamma(n+b+1)*2^(a+b)/((n+a+b+1)*Gamma(n+a+b+1)*(n+1)!*Pn(a,b)'(xi)*P(n+1)(a,b)(xi))
453 
454   Note : the (n+1)/2 first positive points in ascending order only are returned.
455  */
gaussJacobiRule(number_t n,real_t a,real_t b,std::vector<real_t> & points,std::vector<real_t> & weights)456 void gaussJacobiRule(number_t n, real_t a, real_t b, std::vector<real_t>& points, std::vector<real_t>& weights)
457 {
458   if (a==0. && b==0.) { gaussLegendreRule(n, points, weights); }
459   else if (a==2. && b==0.) { gaussJacobi20Rule(n, points, weights); }
460   else error("not_handled");
461 }
462 
463 /*!
464   Compute Gauss-Jacobi n-point formula on [-1, 1]
465   - Points xj are roots of Jacobi Polynomial Pn(2,0) of ordre n (namely orthogonal polynomials of weight (1-x)^2
466   - Weights are - (2n+4)*Gamma(n+3)*Gamma(n+1)*4/((n+3)*Gamma(n+3)*(n+1)!*Pn(2,0)'(xi)*P(n+1)(2,0)(xi))
467 
468   Note : the (n+1)/2 first positive points in ascending order only are returned.
469  */
gaussJacobi20Rule(number_t n,std::vector<real_t> & points,std::vector<real_t> & weights)470 void gaussJacobi20Rule(number_t n, std::vector<real_t>& points, std::vector<real_t>& weights)
471 {
472   std::vector<real_t>::iterator it_p(points.begin()), it_w(weights.begin());
473   switch (n)
474   {
475     case 1:
476       *it_p++ = -0.5;                                *it_w++ = 2.666666666666666666;
477       break;
478     case 2:
479       *it_p++ = -0.754970354689117244;               *it_w++ = 1.860379610028063222;
480       *it_p++ =  0.0883036880224505776;              *it_w++ = 0.8062870566386034446;
481       break;
482     case 3:
483       *it_p++ = -0.854011951853700536;               *it_w++ = 1.257090888519092906;
484       *it_p++ = -0.305992467923296231;               *it_w++ = 1.169970154078928176;
485       *it_p++ =  0.410004419776996766;               *it_w++ = 0.2396056240686455842;
486       break;
487     case 4:
488       *it_p++ = -0.902998901106005341;               *it_w++ = 0.8871073248902238797;
489       *it_p++ = -0.522798524896275390;               *it_w++ = 1.147670318393713663;
490       *it_p++ =  0.0340945902087350046;              *it_w++ = 0.5490710973833846024;
491       *it_p++ =  0.591702835793545726;               *it_w++ = 0.08281792599934452221;
492       break;
493     case 5:
494       *it_p++ = -0.930842120163569817;               *it_w++ = 0.6541182742861673435;
495       *it_p++ = -0.653039358456608554;               *it_w++ = 1.009591695199291904;
496       *it_p++ = -0.220227225868961344;               *it_w++ = 0.7136012897727200017;
497       *it_p++ =  0.268666945261773545;               *it_w++ = 0.2564448057836953539;
498       *it_p++ =  0.702108425894032836;               *it_w++ = 0.03291060162479206367;
499       break;
500     case 6:
501       *it_p++ = -0.948190889812665614;               *it_w++ = 0.50030962181264750327;
502       *it_p++ = -0.736872116684029732;               *it_w++ = 0.85901199789424506110;
503       *it_p++ = -0.395126163954217534;               *it_w++ = 0.75661749398832962811;
504       *it_p++ =  0.0180728263295041679;              *it_w++ = 0.41031656903692968169;
505       *it_p++ =  0.431362254623427838;               *it_w++ = 0.12576237747956041060;
506       *it_p++ =  0.773611232355123733;               *it_w++ = 0.014648606454954381868;
507       break;
508     case 7:
509       *it_p++ = -0.959734452453198986;               *it_w++ = 0.39421201421150496648;
510       *it_p++ = -0.793821941703901970;               *it_w++ = 0.72559059690148915606;
511       *it_p++ = -0.518891747903884927;               *it_w++ = 0.73387042623836203304;
512       *it_p++ = -0.171995710805880507;               *it_w++ = 0.50517102967113038136;
513       *it_p++ =  0.200043026557985860;               *it_w++ = 0.23537769031622891854;
514       *it_p++ =  0.547034493182875002;               *it_w++ = 0.065303405058437556075;
515       *it_p++ =  0.822366333126005527;               *it_w++ = 0.0071415042695136544342;
516       break;
517     case 8:
518       *it_p++ = -0.967804480896157933;               *it_w++ = 0.31823166245352463503;
519       *it_p++ = -0.834198765028697795;               *it_w++ = 0.61454474613778099865;
520       *it_p++ = -0.609049663022520165;               *it_w++ = 0.68227815337550996617;
521       *it_p++ = -0.316696017045595559;               *it_w++ = 0.54757746737322617698;
522       *it_p++ =  0.0111941563689783439;              *it_w++ = 0.32651541110835218546;
523       *it_p++ =  0.339104543648722907;               *it_w++ = 0.13797491024187986251;
524       *it_p++ =  0.631543407166567521;               *it_w++ = 0.035796173704115263926;
525       *it_p++ =  0.857017929919813794;               *it_w++ = 0.0037481422722775780424;
526       break;
527     case 9:
528       *it_p++ = -0.973668228805771019;               *it_w++ = 0.26208116088831777182;
529       *it_p++ = -0.863830940812464825;               *it_w++ = 0.52391629626717305460;
530       *it_p++ = -0.676480966471850716;               *it_w++ = 0.62138855328444403262;
531       *it_p++ = -0.428217823321559205;               *it_w++ = 0.55427516551843767349;
532       *it_p++ = -0.141092709224374415;               *it_w++ = 0.38832502291605206377;
533       *it_p++ =  0.159388112702326253;               *it_w++ = 0.21074624722039868584;
534       *it_p++ =  0.446537143480670867;               *it_w++ = 0.083248932634817896473;
535       *it_p++ =  0.694873684026474640;               *it_w++ = 0.020595189164869784818;
536       *it_p++ =  0.882491728426548423;               *it_w++ = 0.0020900987721557035348;
537       break;
538     case 10:
539       *it_p++ = -0.978063095087651732;               *it_w++ = 0.21947270968077987121;
540       *it_p++ = -0.886203698932684159;               *it_w++ = 0.45018349122246352693;
541       *it_p++ = -0.728099531899542091;               *it_w++ = 0.56055606166933262150;
542       *it_p++ = -0.515437607734952878;               *it_w++ = 0.53961775505150225330;
543       *it_p++ = -0.263984299101324565;               *it_w++ = 0.42307031013571325275;
544       *it_p++ =  0.0076142528297478139;              *it_w++ = 0.27083652013451292698;
545       *it_p++ =  0.279218977309419419;               *it_w++ = 0.13758060037242406444;
546       *it_p++ =  0.530695359096215792;               *it_w++ = 0.051759112548468233764;
547       *it_p++ =  0.743420149148816975;               *it_w++ = 0.012364185557892629036;
548       *it_p++ =  0.901748585281046333;               *it_w++ = 0.0012259202935772867411;
549       break;
550     default:
551       break;
552   }
553 }
554 
555 
556 //! Output of Gauss-Jacobi rule as displayed in function GaussJacobi20Rule
gaussJacobiOutput(number_t nmax,real_t a,real_t b,std::ostream & out)557 void gaussJacobiOutput(number_t nmax, real_t a, real_t b, std::ostream& out)
558 {
559   if (a==0. && b==0.) { gaussLegendreOutput(nmax, out); }
560   else if (a==2. && b==0.) { gaussJacobi20Output(nmax, out); }
561   else error("not_handled");
562 }
563 
564 //! Output of Gauss-Jacobi (2,0) rule as displayed in function GaussJacobi20Rule
gaussJacobi20Output(number_t nmax,std::ostream & out)565 void gaussJacobi20Output(number_t nmax, std::ostream& out)
566 {
567   std::vector<real_t> p(nmax), w(nmax);
568   // Known exactly for n=1, 2, 3 and points only for n=4, 5
569   for (number_t n = 2; n < nmax; n++)
570   {
571     out << "      case " << n << ":" << std::endl;
572     gaussJacobi20Rule(n, p, w);
573     for (number_t i = 0; i < n; i++)
574     {
575       out.setf(std::ios::scientific);
576       out << "         *pi++=" << std::setprecision(19) << p[i] << "L; *wi=" << w[i] << "L;" << std::endl;
577     }
578     out << "         break;" << std::endl;
579   }
580 }
581 
582 } // end of namespace xlifepp
583 
584