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