MGCL V10  V10
MGCL V10
 全て クラス 名前空間 関数 変数 型定義 列挙型 列挙値 フレンド グループ ページ
Gausp.h
1 /********************************************************************/
2 /* Copyright (c) 2015 DG Technologies Inc. and Yuzi Mizuno */
3 /* All rights reserved. */
4 /********************************************************************/
5 #ifndef _MGGausp_HH_
6 #define _MGGausp_HH_
7 
19 template<class func> double mgGausp(
20  func& f,
21  double a,
23  double b,
24  int n=10
25 ){
26  /* Initialized data */
27  const static struct {
28  double e_1;
29  double fill_2[17];
30  double e_3;
31  double fill_4[15];
32  double e_5[2];
33  double fill_6[16];
34  double e_7[2];
35  double fill_8[14];
36  double e_9[3];
37  double fill_10[15];
38  double e_11[3];
39  double fill_12[13];
40  double e_13[4];
41  double fill_14[14];
42  double e_15[4];
43  double fill_16[12];
44  double e_17[5];
45  double fill_18[13];
46  double e_19[5];
47  double fill_20[11];
48  double e_21[6];
49  double fill_22[12];
50  double e_23[6];
51  double fill_24[10];
52  double e_25[7];
53  double fill_26[11];
54  double e_27[7];
55  double fill_28[9];
56  double e_29[8];
57  double fill_30[10];
58  double e_31[8];
59  double fill_32[8];
60  } equiv_7 = { 0., {0}, .5773502691896257, {0}, 0., .7745966692414832,
61  {0}, .3399810435848563, .8611363115940525, {0}, 0.,
62  .5384693101056831, .9061798459386639, {0}, .2386191860831969,
63  .6612093864662644, .932469514203152, {0}, 0.,
64  .4058451513773972, .7415311855993944, .9491079123427584, {0},
65  .1834346424956498, .5255324099163289, .7966664774136267,
66  .9602898564975361, {0}, 0., .3242534234038089,
67  .6133714327005903, .8360311073266357, .968160239507626, {0},
68  .1488743389816312, .4333953941292473, .6794095682990244,
69  .8650633666889845, .9739065285171717, {0}, 0.,
70  .269543155952345, .5190961292068118, .7301520055740493,
71  .8870625997680953, .9782286581460569, {0}, .1252334085114689,
72  .3678314989981802, .5873179542866174, .7699026741943046,
73  .9041172563704747, .9815606342467192, {0}, 0.,
74  .2304583159551348, .4484927510364469, .6423493394403402,
75  .8015780907333098, .9175983992229779, .9841830547185881, {0},
76  .1080549487073436, .3191123689278898, .515248636358154,
77  .6872929048116854, .827201315069765, .9284348836635735,
78  .9862838086968123, {0}, 0., .2011940939974345,
79  .3941513470775634, .5709721726085388, .72441773136017,
80  .8482065834104272, .9372733924007058, .9879925180204854, {0},
81  .09501250983763744, .2816035507792589, .4580167776572274,
82  .6178762444026437, .7554044083550029, .8656312023878317,
83  .9445750230732326, .9894009349916499 };
84 
85 #define c ((double *)&equiv_7)
86 
87  const static struct {
88  double e_1;
89  double fill_2[17];
90  double e_3;
91  double fill_4[15];
92  double e_5[2];
93  double fill_6[16];
94  double e_7[2];
95  double fill_8[14];
96  double e_9[3];
97  double fill_10[15];
98  double e_11[3];
99  double fill_12[13];
100  double e_13[4];
101  double fill_14[14];
102  double e_15[4];
103  double fill_16[12];
104  double e_17[5];
105  double fill_18[13];
106  double e_19[5];
107  double fill_20[11];
108  double e_21[6];
109  double fill_22[12];
110  double e_23[6];
111  double fill_24[10];
112  double e_25[7];
113  double fill_26[11];
114  double e_27[7];
115  double fill_28[9];
116  double e_29[8];
117  double fill_30[10];
118  double e_31[8];
119  double fill_32[8];
120  } equiv_8 = { 2., {0}, .9999999999999998, {0}, .8888888888888888,
121  .5555555555555558, {0}, .6521451548625459, .3478548451374539,
122  {0}, .5688888888888888, .4786286704993666, .2369268850561892,
123  {0}, .4679139345726909, .3607615730481386, .1713244923791703,
124  {0}, .4179591836734694, .3818300505051189, .2797053914892767,
125  .1294849661688699, {0}, .3626837833783622, .3137066458778873,
126  .2223810344533746, .1012285362903763, {0}, .3302393550012598,
127  .3123470770400028, .2606106964029356, .1806481606948573,
128  .08127438836157435, {0}, .2955242247147527, .2692667193099963,
129  .2190863625159819, .1494513491505806, .06667134430868799, {0}
130  , .2729250867779006, .2628045445102467, .2331937645919905,
131  .1862902109277342, .1255803694649046, .05566856711617373, {0},
132  .2491470458134029, .2334925365383548, .2031674267230659,
133  .1600783285433463, .1069393259953185, .04717533638651187, {0},
134  .2325515532308739, .2262831802628972, .2078160475368886,
135  .1781459807619456, .1388735102197873, .09212149983772848,
136  .04048400476531587, {0}, .2152638534631578, .2051984637212956,
137  .1855383974779379, .1572031671581936, .1215185706879031,
138  .08015808715976016, .03511946033175195, {0},
139  .2025782419255613, .1984314853271116, .1861610000155623,
140  .166269205816994, .1395706779261542, .1071592204671719,
141  .07036604748810814, .0307532419961171, {0}, .1894506104550686,
142  .1826034150449236, .1691565193950025, .1495959888165767,
143  .1246289712555339, .09515851168249285, .06225352393864789,
144  .0271524594117541 };
145 
146 #define w ((double *)&equiv_8)
147 
148  if (n < 1 || n > 16) return 0.;
149 
150  double c1 = (b + a) / 2, c2 = (b - a) / 2;
151 
152  int nh;
153  double v;
154  if(n % 2 == 0){
155  nh = n / 2;
156  v = 0.;
157  }else{
158  nh = (n - 1) / 2;
159  v = w[n * 17 - 17] * f(c1);
160  }
161 
162  for (int i = 1; i <= nh; ++i) {
163  double t1 = c1 + c2 * c[i + n * 17 - 17];
164  double t2 = c1 - c2 * c[i + n * 17 - 17];
165  v += w[i + n * 17 - 17] * (f(t1) + f(t2));
166  }
167 
168  v = c2 * v;
169  return v;
170 }
171 
172 #undef w
173 #undef c
174  // end of ALGORITHM group
176 #endif
mgGausp
double mgGausp(func &f, double a, double b, int n=10)
Legendre-Gauss quadratuer formula over (a,b) . The DE formula (double exponential formula) is applied...
Definition: Gausp.h:19