KADATH
leg.cpp
1 /*
2  Copyright 2017 Philippe Grandclement
3 
4  This file is part of Kadath.
5 
6  Kadath is free software: you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation, either version 3 of the License, or
9  (at your option) any later version.
10 
11  Kadath is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with Kadath. If not, see <http://www.gnu.org/licenses/>.
18 */
19 
20 #include <math.h>
21 #include "headcpp.hpp"
22 #include "array.hpp"
23 #define NBR_LEG_MAX 20
24 
25 namespace Kadath {
26 void legendre (int n, double& poly, double& pder, double& polym1, double& pderm1,
27  double& polym2, double& pderm2, double x) {
28 
29 
30  if (n==0) {
31  poly = 1 ;
32  pder = 0 ;
33  }
34  else
35  if (n==1) {
36  polym1 = 1 ;
37  pderm1 = 0 ;
38  poly = x ;
39  pder = 1 ;
40  }
41  else {
42  polym1 = 1 ;
43  pderm1 = 0 ;
44  poly = x ;
45  pder = 1 ;
46  for (int i=1 ; i<n ; i++) {
47  polym2 = polym1 ;
48  pderm2 = pderm1 ;
49  polym1 = poly ;
50  pderm1 = pder ;
51  poly = ((2*i+1)*x*polym1 - i*polym2)/(i+1) ;
52  pder = ((2*i+1)*polym1+(2*i+1)*x*pderm1-i*pderm2)/(i+1) ;
53  }
54  }
55 }
56 
57 double leg (int n, double x) {
58 
59  double p, dp, p1, dp1, p2, dp2 ;
60  legendre (n, p, dp, p1, dp1, p2, dp2, x) ;
61  return p ;
62 }
63 
64 void gauss_lobato_legendre (int n, Array<double>& coloc, Array<double>& weight,
65  double prec = PRECISION, int itemax = 500) {
66  // Verifs :
67  assert (coloc.get_ndim()==1) ;
68  assert (coloc.get_size(0)==n) ;
69  assert (weight.get_ndim()==1) ;
70  assert (weight.get_size(0)==n) ;
71 
72  double x_plus = 1 ;
73  double x_moins = -1 ;
74  double p_plus, dp_plus, p_plus_m1, dp_plus_m1, p_plus_m2, dp_plus_m2 ;
75  double p_moins, dp_moins, p_moins_m1, dp_moins_m1, p_moins_m2, dp_moins_m2 ;
76  double p, dp, p_m1, dp_m1, p_m2, dp_m2 ;
77 
78  legendre (n, p_plus, dp_plus, p_plus_m1, dp_plus_m1, p_plus_m2, dp_plus_m2, x_plus) ;
79  legendre (n, p_moins, dp_moins, p_moins_m1, dp_moins_m1, p_moins_m2, dp_moins_m2, x_moins) ;
80 
81  double det = p_plus_m1*p_moins_m2 - p_moins_m1*p_plus_m2 ;
82  double r_plus = -p_plus ;
83  double r_moins = -p_moins ;
84  double a = (r_plus*p_moins_m2 - r_moins*p_plus_m2)/det ;
85  double b = (r_moins*p_plus_m1 - r_plus*p_moins_m1)/det ;
86 
87  coloc.set(n-1) = 1 ;
88  double dth = M_PI/(2*n+1) ;
89  double cd = cos (2*dth) ;
90  double sd = sin (2*dth) ;
91  double cs = cos(dth) ;
92  double ss = sin(dth) ;
93 
94  int borne_sup = (n%2==0) ?
95  n/2 : (n+1)/2 ;
96 
97  for (int j=1 ; j<borne_sup ; j++) {
98  double x = cs ;
99  bool loop = true ;
100  int ite = 0 ;
101  while (loop) {
102  legendre (n, p, dp, p_m1, dp_m1, p_m2, dp_m2, x) ;
103  double poly = p + a*p_m1 + b*p_m2 ;
104  double pder = dp + a * dp_m1 + b*dp_m2 ;
105  double sum = 0 ;
106  for (int i=0 ; i<j ; i++)
107  sum += 1./(x-coloc(n-i-1)) ;
108 
109  double increm = -poly/(pder-sum*poly) ;
110 
111  x += increm ;
112  ite ++ ;
113  if ((fabs(increm) < prec) || (ite >itemax))
114  loop = false ;
115  }
116  if (ite > itemax) {
117  cout << "Too many iterations..." << endl ;
118  abort() ;
119  }
120  coloc.set(n-j-1) = x ;
121  double auxi = cs*cd-ss*sd ;
122  ss = cs*sd+ss*cd ;
123  cs = auxi ;
124  }
125  if (n%2==1)
126  coloc.set((n-1)/2) = 0 ;
127  // Copy of the symetric ones :
128  for (int i=0 ; i<borne_sup ; i++)
129  coloc.set(i) = - coloc(n-i-1) ;
130 
131  for (int i=0 ; i<n ; i++) {
132  legendre (n-1, p, dp, p_m1, dp_m1, p_m2, dp_m2, coloc(i)) ;
133  weight.set(i) = 2./(n-1)/n/p/p ;
134  }
135 }
136 
137 
138 double get_things_leg (int indice, int n, int thing) {
139 
140  static bool first = true ;
141  static int nbr_done = 0 ;
142  static int n_done[NBR_LEG_MAX] ;
143  static Array<double>** coloc = new Array<double>* [NBR_LEG_MAX] ;
144  static Array<double>** weight = new Array<double>* [NBR_LEG_MAX] ;
145  static Array<double>** gamma = new Array<double>* [NBR_LEG_MAX] ;
146 
147  // Initialisation
148  if (first) {
149  for (int i=0 ; i<NBR_LEG_MAX ; i++) {
150  n_done[i] = 0 ;
151  coloc[i] = 0x0 ;
152  weight[i] = 0x0 ;
153  gamma[i] = 0x0 ;
154  }
155  first = false ;
156  }
157 
158  // Has the computation been done ?
159  int index = -1 ;
160  for (int i=0 ; i<nbr_done ; i++)
161  if (n_done[i] == n)
162  index = i ;
163 
164  if (index == -1) {
165 
166  if (nbr_done>=NBR_LEG_MAX) {
167  cout << "Too many different number of points for Legendre..." << endl ;
168  abort() ;
169  }
170  // Computation must be done :
171  Array<double> col (n) ;
172  Array<double> poids(n) ;
173 
174  gauss_lobato_legendre (n, col, poids) ;
175  coloc[nbr_done] = new Array<double>(col) ;
176  weight[nbr_done] = new Array<double>(poids) ;
177  gamma[nbr_done] = new Array<double>(n) ;
178  (*gamma[nbr_done]) = 0 ;
179  for (int i=0 ; i<n ; i++)
180  for (int j=0 ; j<n ; j++)
181  gamma[nbr_done]->set(i) += pow(leg(i,(*coloc[nbr_done])(j)),2) *
182  (*weight[nbr_done])(j) ;
183  n_done[nbr_done] = n ;
184  nbr_done ++ ;
185  index = nbr_done - 1 ;
186  }
187 
188  switch (thing) {
189  case 0:
190  return (*coloc[index])(indice) ;
191  break ;
192  case 1:
193  return (*weight[index])(indice) ;
194  break ;
195  case 2:
196  return (*gamma[index])(indice) ;
197  break ;
198  default:
199  cout << "Strange error in get_things_legendre..." << endl ;
200  abort() ;
201  }
202 }
203 
204 
205 double coloc_leg (int i, int nbr) {
206  assert ((i>=0) && (i<nbr)) ;
207  return get_things_leg(i, nbr, 0) ;
208 }
209 
210 double weight_leg (int i, int nbr) {
211  assert ((i>=0) && (i<nbr)) ;
212  return get_things_leg(i, nbr, 1) ;
213 }
214 
215 double gamma_leg (int i, int nbr) {
216  assert ((i>=0) && (i<nbr)) ;
217  return get_things_leg(i, nbr, 2) ;
218 }
219 
220  // Case given parity odd or even...
221 
222 double get_things_leg_parity (int indice, int n, int thing) {
223 
224  static bool first = true ;
225  static int nbr_done = 0 ;
226  static int n_done[NBR_LEG_MAX] ;
227  static Array<double>** coloc = new Array<double>* [NBR_LEG_MAX] ;
228  static Array<double>** weight = new Array<double>* [NBR_LEG_MAX] ;
229 
230  // Initialisation
231  if (first) {
232  for (int i=0 ; i<NBR_LEG_MAX ; i++) {
233  n_done[i] = 0 ;
234  coloc[i] = 0x0 ;
235  weight[i] = 0x0 ;
236  }
237  first = false ;
238  }
239 
240  // Has the computation been done ?
241  int index = -1 ;
242  for (int i=0 ; i<nbr_done ; i++)
243  if (n_done[i] == n)
244  index = i ;
245 
246  if (index == -1) {
247 
248  if (nbr_done>=NBR_LEG_MAX) {
249  cout << "Too many number of points for Legendre..." << endl ;
250  abort() ;
251  }
252  // Computation must be done :
253  Array<double> col (2*n-1) ;
254  Array<double> poids(2*n-1) ;
255 
256  gauss_lobato_legendre (2*n-1, col, poids) ;
257  coloc[nbr_done] = new Array<double>(n) ;
258  for (int i=0 ; i<n ; i++)
259  coloc[nbr_done]->set(i) = col(i+n-1) ;
260  weight[nbr_done] = new Array<double>(n) ;
261  weight[nbr_done]->set(0) = poids(n-1)/2. ;
262  for (int i=1 ; i<n ; i++)
263  weight[nbr_done]->set(i) = poids(i+n-1) ;
264  n_done[nbr_done] = n ;
265  nbr_done ++ ;
266  index = nbr_done - 1 ;
267  }
268 
269  switch (thing) {
270  case 0:
271  return (*coloc[index])(indice) ;
272  break ;
273  case 1:
274  return (*weight[index])(indice) ;
275  break ;
276  default:
277  cout << "Strange error in get_things_legendre..." << endl ;
278  abort() ;
279  }
280 }
281 
282 
283 double coloc_leg_parity (int i, int nbr) {
284  assert ((i>=0) && (i<nbr)) ;
285  return get_things_leg_parity(i, nbr, 0) ;
286 }
287 
288 double weight_leg_parity (int i, int nbr) {
289  assert ((i>=0) && (i<nbr)) ;
290  return get_things_leg_parity(i, nbr, 1) ;
291 }
292 
293 
294 double gamma_leg_even (int nn, int nbr) {
295  assert ((nn>=0) && (nn<nbr)) ;
296  static bool first = true ;
297  static int nbr_done = 0 ;
298  static int n_done[NBR_LEG_MAX] ;
299  static Array<double>** gamma = new Array<double>* [NBR_LEG_MAX] ;
300 
301  if (first) {
302  for (int i=0 ; i<NBR_LEG_MAX ; i++) {
303  n_done[i] = 0 ;
304  gamma[i] = 0x0 ;
305  }
306  first = false ;
307  }
308 
309  // Has the computation been done ?
310  int index = -1 ;
311  for (int i=0 ; i<nbr_done ; i++)
312  if (n_done[i] == nbr)
313  index = i ;
314 
315  if (index == -1) {
316 
317  if (nbr_done>=NBR_LEG_MAX) {
318  cout << "Too many number of points for Legendre..." << endl ;
319  abort() ;
320  }
321  gamma[nbr_done] = new Array<double> (nbr) ;
322  *gamma[nbr_done] = 0 ;
323  for (int i=0 ; i<nbr ; i++)
324  for (int j=0 ; j<nbr ; j++)
325  gamma[nbr_done]->set(i) += pow(leg(i*2,coloc_leg_parity(j,nbr)),2) * weight_leg_parity(j, nbr) ;
326  n_done[nbr_done] = nbr ;
327  nbr_done ++ ;
328  index = nbr_done - 1 ;
329  }
330 
331  return (*gamma[index])(nn) ;
332 }
333 
334 
335 double gamma_leg_odd (int nn, int nbr) {
336  assert ((nn>=0) && (nn<nbr)) ;
337  static bool first = true ;
338  static int nbr_done = 0 ;
339  static int n_done[NBR_LEG_MAX] ;
340  static Array<double>** gamma = new Array<double>* [NBR_LEG_MAX] ;
341 
342  if (first) {
343  for (int i=0 ; i<NBR_LEG_MAX ; i++) {
344  n_done[i] = 0 ;
345  gamma[i] = 0x0 ;
346  }
347  first = false ;
348  }
349 
350  // Has the computation been done ?
351  int index = -1 ;
352  for (int i=0 ; i<nbr_done ; i++)
353  if (n_done[i] == nbr)
354  index = i ;
355 
356  if (index == -1) {
357 
358  if (nbr_done>=NBR_LEG_MAX) {
359  cout << "Too many number of points for Legendre..." << endl ;
360  abort() ;
361  }
362  gamma[nbr_done] = new Array<double> (nbr) ;
363  *gamma[nbr_done] = 0 ;
364  for (int i=0 ; i<nbr ; i++)
365  for (int j=0 ; j<nbr ; j++)
366  gamma[nbr_done]->set(i) += pow(leg(i*2+1,coloc_leg_parity(j,nbr)),2) * weight_leg_parity(j, nbr) ;
367  n_done[nbr_done] = nbr ;
368  nbr_done ++ ;
369  index = nbr_done - 1 ;
370  }
371 
372  return (*gamma[index])(nn) ;
373 }
374 }
375