% This example code does adaptive integration of single variable functions. % % ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; % % TRAPEZOIDAL INTEGRATOR % % ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; % % Used as a step in the adaptive method % function trapezoidal(func,a,b,n) = let delta = (b-a)/float(n) in delta * sum({func(delta * float(i) + a) : i in index(n)}) $ % ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; % % SIMPLE ADAPTIVE INTEGRATOR % % ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; % % Some constants % accuracy1 = 1e-3; % Error as a fraction of the value of the integral (approximate) % points1 = 1000; % Number of points used in trapezoidal integral % % Simple adaptive integrator. func -- a function, must have type (float -> float) (a,b) -- an interval, must both be floats % function adaptive1(func,(a,b)) = let integral = trapezoidal(func,a,b,points1); integral2 = trapezoidal(func,a,b,points1*2); mid = (a+b)/2.0 in if abs(integral2 - integral) > abs(accuracy1*integral) then % Parallel recursive calls on the two halves of the interval % sum({adaptive1(func,interval): interval in [(a,mid),(mid,b)]}) else integral2 $ function integrate1(func,a,b) = adaptive1(func,a,b); % ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; % % IMPROVED ADAPTIVE INTEGRATOR % % ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; % % Constants % accuracy2 = 1e-4; points2 = 1000; % Number of points used in trapezoidal integral % conv_depth = 7; % Depth of recursion at which start making serial calls % max_depth = 9; % Depth of recursion at which it quits % % This version has two improvements over the first version. 1) It has a bound on how deep the recursion can go. This can prevent infinite loops at singularities. 2) It starts making serial calls when recursion reaches a certain depth. This reduces parallelism and can prevent memory overflows. % function adaptive2(func,depth,a,b) = let integral = trapezoidal(func,a,b,points2); integral2 = trapezoidal(func,a,b,points2*2); mid = (a+b)/2.0 in if (abs(integral2 - integral) > abs(accuracy2*integral)) and (depth < max_depth) then if depth < conv_depth then sum({adaptive2(func,depth+1,interval): interval in [(a,mid),(mid,b)]}) else adaptive2(func,depth+1,a,mid) + adaptive2(func,depth+1,mid,b) else integral2 $ function integrate2(func,a,b) = adaptive2(func,0,a,b);