_b_e_g_i_n _c_o_m_m_e_n_t HAVIE INTEGRATOR.
ALGORITHM 257, Robert N. Kubik,
CACM 8 (1965) 381;
_r_e_a_l a,b,eps,mask,y,answer; _i_n_t_e_g_e_r count;
_r_e_a_l _p_r_o_c_e_d_u_r_e havieintegrator(x,a,b,eps,integrand,m);
_v_a_l_u_e a,b,eps,m;
_i_n_t_e_g_e_r m; _r_e_a_l integrand,x,a,b,eps;
_c_o_m_m_e_n_t This algorithm performs numerical integration of
definite integrals using an equidistant sampling of the
function and repeated halving of the sampling interval.
Each halving allows the calculation of a trapezium and
a tangent formula on a finer grid, but also the calcul-
ation of several higher order formulas which are defined
implicitly. The two families of approximate solutions
will normally bracket the value of the integral and from
these convergence is tested on each of the several orders
of approximation. The algorithm is based on a private
communication from F. Haavie of the Institutt for Atom-
energi Kjeller Research Establishment, Norway. A Fortran
version is in use on the Philco-2000. ...;
_b_e_g_i_n _r_e_a_l h,endpts,sumt,sumu,d;
_i_n_t_e_g_e_r i,j,k,n;
_r_e_a_l _a_r_r_a_y t,u,tprev,uprev[1:m];
x:= a; endpts:= integrand; count:= 1;
x:= b; endpts:= 0.5 * (integrand + endpts); count:= count + 1;
sumt:= 0.0; i:= n:= 1; h:= b - a;
estimate:
t[1]:= h * (endpts + sumt); sumu:= 0.0;
_c_o_m_m_e_n_t t[1] = h*(0.5*f[0]+f[1]+f[2]+...+0.5*f[2^(i-1)]);
x:= a - h/2.0;
_f_o_r j:= 1 _s_t_e_p 1 _u_n_t_i_l n _d_o
_b_e_g_i_n x:= x + h; sumu:= sumu + integrand; count:= count + 1 _e_n_d;
u[1]:= h * sumu; k:= 1;
_c_o_m_m_e_n_t u[1] = h*(f[1/2]+f[3/2]+...f[(2^i-1)/2],
k corresponds to approximate solution with truncation
error term of order 2k;
test:
_i_f abs(t[k]-u[k]) _< eps
_t_h_e_n _b_e_g_i_n havieintegrator:= 0.5 * (t[k] + u[k]);
_g_o_t_o exit
_e_n_d;
_i_f k |= i
_t_h_e_n _b_e_g_i_n d:= 2 |^ (2*k);
t[k+1]:= (d * t[k] - tprev[k]) / (d - 1.0);
tprev[k]:= t[k];
u[k+1]:= (d * u[k] - uprev[k]) / (d - 1.0);
uprev[k]:= u[k];
_c_o_m_m_e_n_t This implicit formulation of the higher
order integration formulas is given in
[ROMBERG, W. ...;
k:= k + 1;
_i_f k = m
_t_h_e_n _b_e_g_i_n havieintegrator:= mask;
_g_o_t_o exit
_e_n_d;
_g_o_t_o test
_e_n_d;
h:= h / 2.0; sumt:= sumt + sumu;
tprev[k]:= t[k]; uprev[k]:= u[k];
i:= i + 1; n:= 2 * n;
_g_o_t_o estimate;
exit: NLCR; NLCR;
PRINTTEXT(|<i: |>); ABSFIXT(3,0,i);
PRINTTEXT(|<k: |>); ABSFIXT(3,0,k);
PRINTTEXT(|< count:|>); ABSFIXT(7,0,count)
_e_n_d havieintegrator;
_c_o_m_m_e_n_t Following is a driver program to test havieintegrator;
a:= 0.0; b:= 1.5707963; eps:= 0.00001; mask:= 9.99;
answer:= havieintegrator(y,a,b,eps,cos(y),12);
NLCR; PRINTTEXT(|<integral(0,pi/2,cos(x)) = |>);
FLOT(10,2,answer);
a:= 0.0; b:= 4.3;
answer:= havieintegrator(y,a,b,eps,exp(-y*y),12);
NLCR; PRINTTEXT(|<integral(0,4.3,exp(-x*x)) = |>);
FLOT(10,2,answer);
a:= 1.0; b:= 10.0;
answer:= havieintegrator(y,a,b,eps,ln(y),12);
NLCR; PRINTTEXT(|<integral(1,10,ln(x)) = |>);
FLOT(10,2,answer);
a:= 0.0; b:= 20.0;
answer:= havieintegrator(y,a,b,eps,y|^(1/2)/(exp(y-4)+1),20);
NLCR; PRINTTEXT(|<integral(1,20,x^(1/2)/(exp(x-4)+1)) = |>);
FLOT(10,2,answer);
_e_n_d