1 /* integration/qmomof.c
2 *
3 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
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 (at
8 * your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful, but
11 * WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
18 */
19
20 #include "gsl__config.h"
21 #include <stdlib.h>
22 #include "gsl_integration.h"
23 #include "gsl_errno.h"
24
25 static void
26 compute_moments (double par, double * cheb);
27
28 static int
29 dgtsl (size_t n, double *c, double *d, double *e, double *b);
30
31 gsl_integration_qawo_table *
gsl_integration_qawo_table_alloc(double omega,double L,enum gsl_integration_qawo_enum sine,size_t n)32 gsl_integration_qawo_table_alloc (double omega, double L,
33 enum gsl_integration_qawo_enum sine,
34 size_t n)
35 {
36 gsl_integration_qawo_table *t;
37 double * chebmo;
38
39 if (n == 0)
40 {
41 GSL_ERROR_VAL ("table length n must be positive integer",
42 GSL_EDOM, 0);
43 }
44
45 t = (gsl_integration_qawo_table *)
46 malloc (sizeof (gsl_integration_qawo_table));
47
48 if (t == 0)
49 {
50 GSL_ERROR_VAL ("failed to allocate space for qawo_table struct",
51 GSL_ENOMEM, 0);
52 }
53
54 chebmo = (double *) malloc (25 * n * sizeof (double));
55
56 if (chebmo == 0)
57 {
58 free (t);
59 GSL_ERROR_VAL ("failed to allocate space for chebmo block",
60 GSL_ENOMEM, 0);
61 }
62
63 t->n = n;
64 t->sine = sine;
65 t->omega = omega;
66 t->L = L;
67 t->par = 0.5 * omega * L;
68 t->chebmo = chebmo;
69
70 /* precompute the moments */
71
72 {
73 size_t i;
74 double scale = 1.0;
75
76 for (i = 0 ; i < t->n; i++)
77 {
78 compute_moments (t->par * scale, t->chebmo + 25*i);
79 scale *= 0.5;
80 }
81 }
82
83 return t;
84 }
85
86 int
gsl_integration_qawo_table_set(gsl_integration_qawo_table * t,double omega,double L,enum gsl_integration_qawo_enum sine)87 gsl_integration_qawo_table_set (gsl_integration_qawo_table * t,
88 double omega, double L,
89 enum gsl_integration_qawo_enum sine)
90 {
91 t->omega = omega;
92 t->sine = sine;
93 t->L = L;
94 t->par = 0.5 * omega * L;
95
96 /* recompute the moments */
97
98 {
99 size_t i;
100 double scale = 1.0;
101
102 for (i = 0 ; i < t->n; i++)
103 {
104 compute_moments (t->par * scale, t->chebmo + 25*i);
105 scale *= 0.5;
106 }
107 }
108
109 return GSL_SUCCESS;
110 }
111
112
113 int
gsl_integration_qawo_table_set_length(gsl_integration_qawo_table * t,double L)114 gsl_integration_qawo_table_set_length (gsl_integration_qawo_table * t,
115 double L)
116 {
117 /* return immediately if the length is the same as the old length */
118
119 if (L == t->L)
120 return GSL_SUCCESS;
121
122 /* otherwise reset the table and compute the new parameters */
123
124 t->L = L;
125 t->par = 0.5 * t->omega * L;
126
127 /* recompute the moments */
128
129 {
130 size_t i;
131 double scale = 1.0;
132
133 for (i = 0 ; i < t->n; i++)
134 {
135 compute_moments (t->par * scale, t->chebmo + 25*i);
136 scale *= 0.5;
137 }
138 }
139
140 return GSL_SUCCESS;
141 }
142
143
144 void
gsl_integration_qawo_table_free(gsl_integration_qawo_table * t)145 gsl_integration_qawo_table_free (gsl_integration_qawo_table * t)
146 {
147 free (t->chebmo);
148 free (t);
149 }
150
151 static void
compute_moments(double par,double * chebmo)152 compute_moments (double par, double *chebmo)
153 {
154 double v[28], d[25], d1[25], d2[25];
155
156 const size_t noeq = 25;
157
158 const double par2 = par * par;
159 const double par4 = par2 * par2;
160 const double par22 = par2 + 2.0;
161
162 const double sinpar = sin (par);
163 const double cospar = cos (par);
164
165 size_t i;
166
167 /* compute the chebyschev moments with respect to cosine */
168
169 double ac = 8 * cospar;
170 double as = 24 * par * sinpar;
171
172 v[0] = 2 * sinpar / par;
173 v[1] = (8 * cospar + (2 * par2 - 8) * sinpar / par) / par2;
174 v[2] = (32 * (par2 - 12) * cospar
175 + (2 * ((par2 - 80) * par2 + 192) * sinpar) / par) / par4;
176
177 if (fabs (par) <= 24)
178 {
179 /* compute the moments as the solution of a boundary value
180 problem using the asyptotic expansion as an endpoint */
181
182 double an2, ass, asap;
183 double an = 6;
184 size_t k;
185
186 for (k = 0; k < noeq - 1; k++)
187 {
188 an2 = an * an;
189 d[k] = -2 * (an2 - 4) * (par22 - 2 * an2);
190 d2[k] = (an - 1) * (an - 2) * par2;
191 d1[k + 1] = (an + 3) * (an + 4) * par2;
192 v[k + 3] = as - (an2 - 4) * ac;
193 an = an + 2.0;
194 }
195
196 an2 = an * an;
197
198 d[noeq - 1] = -2 * (an2 - 4) * (par22 - 2 * an2);
199 v[noeq + 2] = as - (an2 - 4) * ac;
200 v[3] = v[3] - 56 * par2 * v[2];
201
202 ass = par * sinpar;
203 asap = (((((210 * par2 - 1) * cospar - (105 * par2 - 63) * ass) / an2
204 - (1 - 15 * par2) * cospar + 15 * ass) / an2
205 - cospar + 3 * ass) / an2
206 - cospar) / an2;
207 v[noeq + 2] = v[noeq + 2] - 2 * asap * par2 * (an - 1) * (an - 2);
208
209 dgtsl (noeq, d1, d, d2, v + 3);
210
211 }
212 else
213 {
214 /* compute the moments by forward recursion */
215 size_t k;
216 double an = 4;
217
218 for (k = 3; k < 13; k++)
219 {
220 double an2 = an * an;
221 v[k] = ((an2 - 4) * (2 * (par22 - 2 * an2) * v[k - 1] - ac)
222 + as - par2 * (an + 1) * (an + 2) * v[k - 2])
223 / (par2 * (an - 1) * (an - 2));
224 an = an + 2.0;
225 }
226 }
227
228
229 for (i = 0; i < 13; i++)
230 {
231 chebmo[2 * i] = v[i];
232 }
233
234 /* compute the chebyschev moments with respect to sine */
235
236 v[0] = 2 * (sinpar - par * cospar) / par2;
237 v[1] = (18 - 48 / par2) * sinpar / par2 + (-2 + 48 / par2) * cospar / par;
238
239 ac = -24 * par * cospar;
240 as = -8 * sinpar;
241
242 if (fabs (par) <= 24)
243 {
244 /* compute the moments as the solution of a boundary value
245 problem using the asyptotic expansion as an endpoint */
246
247 size_t k;
248 double an2, ass, asap;
249 double an = 5;
250
251 for (k = 0; k < noeq - 1; k++)
252 {
253 an2 = an * an;
254 d[k] = -2 * (an2 - 4) * (par22 - 2 * an2);
255 d2[k] = (an - 1) * (an - 2) * par2;
256 d1[k + 1] = (an + 3) * (an + 4) * par2;
257 v[k + 2] = ac + (an2 - 4) * as;
258 an = an + 2.0;
259 }
260
261 an2 = an * an;
262
263 d[noeq - 1] = -2 * (an2 - 4) * (par22 - 2 * an2);
264 v[noeq + 1] = ac + (an2 - 4) * as;
265 v[2] = v[2] - 42 * par2 * v[1];
266
267 ass = par * cospar;
268 asap = (((((105 * par2 - 63) * ass - (210 * par2 - 1) * sinpar) / an2
269 + (15 * par2 - 1) * sinpar
270 - 15 * ass) / an2 - sinpar - 3 * ass) / an2 - sinpar) / an2;
271 v[noeq + 1] = v[noeq + 1] - 2 * asap * par2 * (an - 1) * (an - 2);
272
273 dgtsl (noeq, d1, d, d2, v + 2);
274
275 }
276 else
277 {
278 /* compute the moments by forward recursion */
279 size_t k;
280 double an = 3;
281 for (k = 2; k < 12; k++)
282 {
283 double an2 = an * an;
284 v[k] = ((an2 - 4) * (2 * (par22 - 2 * an2) * v[k - 1] + as)
285 + ac - par2 * (an + 1) * (an + 2) * v[k - 2])
286 / (par2 * (an - 1) * (an - 2));
287 an = an + 2.0;
288 }
289 }
290
291 for (i = 0; i < 12; i++)
292 {
293 chebmo[2 * i + 1] = v[i];
294 }
295
296 }
297
298 static int
dgtsl(size_t n,double * c,double * d,double * e,double * b)299 dgtsl (size_t n, double *c, double *d, double *e, double *b)
300 {
301 /* solves a tridiagonal matrix A x = b
302
303 c[1 .. n - 1] subdiagonal of the matrix A
304 d[0 .. n - 1] diagonal of the matrix A
305 e[0 .. n - 2] superdiagonal of the matrix A
306
307 b[0 .. n - 1] right hand side, replaced by the solution vector x */
308
309 size_t k;
310
311 c[0] = d[0];
312
313 if (n == 0)
314 {
315 return GSL_SUCCESS;
316 }
317
318 if (n == 1)
319 {
320 b[0] = b[0] / d[0] ;
321 return GSL_SUCCESS;
322 }
323
324 d[0] = e[0];
325 e[0] = 0;
326 e[n - 1] = 0;
327
328 for (k = 0; k < n - 1; k++)
329 {
330 size_t k1 = k + 1;
331
332 if (fabs (c[k1]) >= fabs (c[k]))
333 {
334 {
335 double t = c[k1];
336 c[k1] = c[k];
337 c[k] = t;
338 };
339 {
340 double t = d[k1];
341 d[k1] = d[k];
342 d[k] = t;
343 };
344 {
345 double t = e[k1];
346 e[k1] = e[k];
347 e[k] = t;
348 };
349 {
350 double t = b[k1];
351 b[k1] = b[k];
352 b[k] = t;
353 };
354 }
355
356 if (c[k] == 0)
357 {
358 return GSL_FAILURE ;
359 }
360
361 {
362 double t = -c[k1] / c[k];
363
364 c[k1] = d[k1] + t * d[k];
365 d[k1] = e[k1] + t * e[k];
366 e[k1] = 0;
367 b[k1] = b[k1] + t * b[k];
368 }
369
370 }
371
372 if (c[n - 1] == 0)
373 {
374 return GSL_FAILURE;
375 }
376
377
378 b[n - 1] = b[n - 1] / c[n - 1];
379
380 b[n - 2] = (b[n - 2] - d[n - 2] * b[n - 1]) / c[n - 2];
381
382 for (k = n ; k > 2; k--)
383 {
384 size_t kb = k - 3;
385 b[kb] = (b[kb] - d[kb] * b[kb + 1] - e[kb] * b[kb + 2]) / c[kb];
386 }
387
388 return GSL_SUCCESS;
389 }
390