-
Notifications
You must be signed in to change notification settings - Fork 105
/
Copy pathadapt.h
executable file
·86 lines (83 loc) · 2.43 KB
/
adapt.h
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
struct Adapt {
Doub TOL,toler;
static const Doub alpha,beta,x1,x2,x3,x[12];
bool terminate,out_of_tolerance;
Adapt(Doub tol);
template <class T>
Doub integrate(T &func, const Doub a, const Doub b);
template <class T>
Doub adaptlob(T &func, const Doub a, const Doub b, const Doub fa,
const Doub fb, const Doub is);
};
Adapt::Adapt(Doub tol) : TOL(tol),terminate(true),out_of_tolerance(false)
{
const Doub EPS=numeric_limits<Doub>::epsilon();
if (TOL < 10.0*EPS)
TOL=10.0*EPS;
}
template <class T>
Doub Adapt::integrate(T &func, const Doub a, const Doub b)
{
Doub m,h,fa,fb,i1,i2,is,erri1,erri2,r,y[13];
m=0.5*(a+b);
h=0.5*(b-a);
fa=y[0]=func(a);
fb=y[12]=func(b);
for (Int i=1;i<12;i++)
y[i]=func(m+x[i]*h);
i2=(h/6.0)*(y[0]+y[12]+5.0*(y[4]+y[8]));
i1=(h/1470.0)*(77.0*(y[0]+y[12])+432.0*(y[2]+y[10])+
625.0*(y[4]+y[8])+672.0*y[6]);
is=h*(0.0158271919734802*(y[0]+y[12])+0.0942738402188500*
(y[1]+y[11])+0.155071987336585*(y[2]+y[10])+
0.188821573960182*(y[3]+y[9])+0.199773405226859*
(y[4]+y[8])+0.224926465333340*(y[5]+y[7])+
0.242611071901408*y[6]);
erri1=abs(i1-is);
erri2=abs(i2-is);
r=(erri2 != 0.0) ? erri1/erri2 : 1.0;
toler=(r > 0.0 && r < 1.0) ? TOL/r : TOL;
if (is == 0.0)
is=b-a;
is=abs(is);
return adaptlob(func,a,b,fa,fb,is);
}
template <class T>
Doub Adapt::adaptlob(T &func, const Doub a, const Doub b, const Doub fa,
const Doub fb, const Doub is)
{
Doub m,h,mll,ml,mr,mrr,fmll,fml,fm,fmrr,fmr,i1,i2;
m=0.5*(a+b);
h=0.5*(b-a);
mll=m-alpha*h;
ml=m-beta*h;
mr=m+beta*h;
mrr=m+alpha*h;
fmll=func(mll);
fml=func(ml);
fm=func(m);
fmr=func(mr);
fmrr=func(mrr);
i2=h/6.0*(fa+fb+5.0*(fml+fmr));
i1=h/1470.0*(77.0*(fa+fb)+432.0*(fmll+fmrr)+625.0*(fml+fmr)+672.0*fm);
if (abs(i1-i2) <= toler*is || mll <= a || b <= mrr) {
if ((mll <= a || b <= mrr) && terminate) {
out_of_tolerance=true;
terminate=false;
}
return i1;
}
else
return adaptlob(func,a,mll,fa,fmll,is)+
adaptlob(func,mll,ml,fmll,fml,is)+
adaptlob(func,ml,m,fml,fm,is)+
adaptlob(func,m,mr,fm,fmr,is)+
adaptlob(func,mr,mrr,fmr,fmrr,is)+
adaptlob(func,mrr,b,fmrr,fb,is);
}
const Doub Adapt::alpha=sqrt(2.0/3.0);
const Doub Adapt::beta=1.0/sqrt(5.0);
const Doub Adapt::x1=0.942882415695480;
const Doub Adapt::x2=0.641853342345781;
const Doub Adapt::x3=0.236383199662150;
const Doub Adapt::x[12]={0,-x1,-alpha,-x2,-beta,-x3,0.0,x3,beta,x2,alpha,x1};