/* deint.c -- integration of function by double exponential transform * from Applied Numerical Methods with Software, by Nakamura. * * (C) 2001, C. Bond. All rights reserved. */ #include #include #include double deint(double (*f)(double),int n,double a,double b) { double h,ss,pi,z,exz,hcos,hsin,s,dxdz,x,w; int k; n = n/2; h = 5.0/n; /* 5.0 is rough limit of K in exp(exp(K)) */ pi = M_PI_4; ss = 0.0; for (k=-n;k<=n;k++) { z = h * (double)k; exz = exp(z); hcos = exz+1.0/exz; hsin = exz-1.0/exz; s = exp(pi*hsin); w = s + 1.0/s; x = (b*s+a/s)/w; /* transformed abscissa */ if (x != a && x != b) { dxdz = 2*(b-a)*pi*hcos/(w*w); /* transformed weight */ ss += h * f(x)*dxdz; } } return ss; } double deint2(double (*f)(double),int n,double a,double b) { double h,ss,pi,z,exz,hcos,hsin,s,dxdz,x1,x2,w; int k; n = n/2; h = 5.0/n; /* 5.0 is rough limit of K in exp(exp(K)) */ pi = M_PI_4; ss = 0.5*f((a+b)/2); for (k=-n;k<0;k++) { z = h*(double)k; exz = exp(z); hcos = exz+1.0/exz; hsin = exz-1.0/exz; s = exp(pi*hsin); w = s + 1.0/s; dxdz = hcos/(w*w); /* transformed weight */ x1 = (b*s+a/s)/w; /* transformed abscissa */ x2 = (a*s+b/s)/w; /* " " */ if (x1 != a && x1 != b) ss += dxdz * f(x1); if (x2 != a && x2 != b) ss += dxdz * f(x2); } return 2*(b-a)*pi*h*ss; } double deint3(double (*f)(double),int n,double a,double b) { double h,ss,pi,z,exz,hcos,hsin,s,dxdz,x,w; int k; n = n/2; h = 5.0/n; /* 5.0 is rough limit of K in exp(exp(K)) */ pi = M_PI_4; ss = 0.0; for (k=-n;k<=n;k++) { z = h * (double)k; exz = exp(z); hcos = exz+1.0/exz; hsin = exz-1.0/exz; s = exp(pi*hsin); w = s + 1.0/s; x = (b*s+a/s)/w; /* transformed abscissa */ if (x != a && x != b) { dxdz = hcos/(w*w); /* transformed weight */ ss += f(x) * dxdz; } } printf ("1st subtotal: %15.15lf\n", h*2*(b-a)*pi*ss); for (k=-n;k<=n;k++) { z = h * ((double)(k)+0.5); /* use in-between points this pass */ exz = exp(z); hcos = exz+1.0/exz; hsin = exz-1.0/exz; s = exp(pi*hsin); w = s + 1.0/s; x = (b*s+a/s)/w; /* transformed abscissa */ if (x != a && x != b) { dxdz = hcos/(w*w); /* transformed weight */ ss += f(x) * dxdz; } } printf ("2nd subtotal: %15.15lf\n", h*(b-a)*pi*ss); return h*(b-a)*pi*ss; } /* Progressive version */ double deint4(double (*f)(double),int n,double a,double b) { double h,ss,ss1,ss2,pi,z,exz,hcos,hsin,s,dxdz,x1,x2,w; int k,l; n = n/2; h = 5.0/n; /* 5.0 is rough limit of K in exp(exp(K)) */ pi = M_PI_4; ss = 0.5*f((a+b)/2); ss1 = 0.0; ss2 = 0.0; for (k=-n;k<0;k++) { z = h*(double)k; exz = exp(z); hcos = exz+1.0/exz; hsin = exz-1.0/exz; s = exp(pi*hsin); w = s + 1.0/s; dxdz = hcos/(w*w); /* transformed weight */ x1 = (b*s+a/s)/w; /* transformed abscissa */ x2 = (a*s+b/s)/w; /* " " */ if (x1 != a && x1 != b) ss1 += dxdz * f(x1); if (x2 != a && x2 != b) ss2 += dxdz * f(x2); } printf("1, %15.15lf, %15.15lf\n",2*(b-a)*pi*h*ss1,2*(b-a)*pi*h*ss2); for (l=0;l<3;l++) { for (k=-n;k<0;k++) { z = h*((double)k+0.5); exz = exp(z); hcos = exz+1.0/exz; hsin = exz-1.0/exz; s = exp(pi*hsin); w = s + 1.0/s; dxdz = hcos/(w*w); /* transformed weight */ x1 = (b*s+a/s)/w; /* transformed abscissa */ x2 = (a*s+b/s)/w; /* " " */ if (x1 != a && x1 != b) ss1 += dxdz * f(x1); if (x2 != a && x2 != b) ss2 += dxdz * f(x2); } printf("%d, %15.15lf, %15.15lf\n",l+2,(b-a)*pi*h*ss1,(b-a)*pi*h*ss2); n *= 2; h /= 2; } return 2*(b-a)*pi*h*(ss+ss1+ss2); }