-
Notifications
You must be signed in to change notification settings - Fork 105
/
dfridr.h
executable file
·31 lines (31 loc) · 781 Bytes
/
dfridr.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
template<class T>
Doub dfridr(T &func, const Doub x, const Doub h, Doub &err)
{
const Int ntab=10;
const Doub con=1.4, con2=(con*con);
const Doub big=numeric_limits<Doub>::max();
const Doub safe=2.0;
Int i,j;
Doub errt,fac,hh,ans;
MatDoub a(ntab,ntab);
if (h == 0.0) throw("h must be nonzero in dfridr.");
hh=h;
a[0][0]=(func(x+hh)-func(x-hh))/(2.0*hh);
err=big;
for (i=1;i<ntab;i++) {
hh /= con;
a[0][i]=(func(x+hh)-func(x-hh))/(2.0*hh);
fac=con2;
for (j=1;j<=i;j++) {
a[j][i]=(a[j-1][i]*fac-a[j-1][i-1])/(fac-1.0);
fac=con2*fac;
errt=MAX(abs(a[j][i]-a[j-1][i]),abs(a[j][i]-a[j-1][i-1]));
if (errt <= err) {
err=errt;
ans=a[j][i];
}
}
if (abs(a[i][i]-a[i-1][i-1]) >= safe*err) break;
}
return ans;
}