Description of Contents of cinte.tgz
( Self-Validating Integration and Approximation
of Piecewise Analytic Functions )
The package is programmed in PASCAL-XSC.
You can compile and use it (at least) with a GNU C-compiler
with library version 2.7.2.1.
The package consists of the files
- makefile :
updates the modules
if necessary
- g_form.p, formeln.p :
nodes and coefficients of Gaussian quadrature formulas
- arith.p :
complex interval arithmetic for functions that are real-valued
for real arguments.
- bi_ari.p : "bounded interval"-arithmetic;
an arithmetic for the calculation of bounds for a function. Also
checks boundedness by a machine number before application of the function,
i.e., instead of exiting due to over- or underflow.
- ade_ari.p : "analytic, differentiable,
evaluation"-arithmetic; an arithmetic that calculates
enclosures for
the
function (real), the first derivative (real) and the
complex function values for interval data. It checks
the boundedness, differentiability and analyticity before
applying the function.
- ade_comf.p :
combination of arithmetics in
i_ari, bi_ari
and ade_ari in one data structure (slower for simple real interval
evaluations)
- cinte.p :
"integration using complex
interval arithmetic"; contains four integration routines
that differ in the parameter list.
- cinte_comf.p :
"comfortable
integration using complex
interval arithmetic"; contains four
integration routines
that differ in the parameter list. Is slower but more comfortable than
cinte.p
- test1.p : test example with several integrands.
Some of them are
also contained in the corresponding preprint (cf.
Software-home-page)
- test2.p : test examples from
QUADPACK
(cf. corresponding preprint from my
Software-home-page)
- test3.p : equivalent to test1.p but using
cinte_comf.p
Use of cinte.p:
Available integration routines in cinte.p are
GLOBAL FUNCTION INTEGRATE(function f(x: interval): interval;
function fb(xx : interval):b_interval;
function fa(xx : interval; r:real):analyticd;
a,b, eps:real): interval;
GLOBAL FUNCTION INTEGRATE(function f(x: interval): interval;
function fb(xx : interval):b_interval;
function fa(xx : interval; r:real):analyticd;
a,b, eps,except:real): interval;
GLOBAL FUNCTION INTEGRATE(function f(x: interval): interval;
function fb(xx : interval):b_interval;
function fa(xx : interval; r:real):analyticd;
a,b, eps; n:integer): interval;
GLOBAL FUNCTION INTEGRATE(function f(x: interval): interval;
function fb(xx : interval):b_interval;
function fa(xx : interval; r:real):analyticd;
a,b, eps,except:real; n:integer): interval;
where
- f,fb,fa : routines representing the same function (see Example
below)
- a,b : lower and upper integration bound
- eps : required precision. This is usually an absolute
precision. If function values become too large, however, it is switched
to a relative precision.
- n : number of used function values
- except : bounds for widths of intervals,
where f is not investigated
further with respect to boundedness if boundedness is not already proved.
The use of the version with parameter except is
meant for possibly unbounded integrands (try
f2 in test1.p or test3.p )
As an example
and in order to understand the meaning of the parameters,
consider the relevant program text from test1.p
program test (input, output);
use cinte;
Now, you can use the integration routines
....
function fa1(xx : interval; r:real):analyticd;
var x : analyticd;
begin
x := construct(xx, r);
fa1:= sin(sqr(x))*exp(x);
end;
function fb1(xx : interval):b_interval;
var x : b_interval;
begin
x := construct(xx);
fb1:= sin(sqr(x))*exp(x);
end;
function f1(x : interval): interval;
begin
f1 := sin(sqr(x))*exp(x);
end;
....
In this faster version, the user is required to provide functions for
different types of arguments.
- The first function fa1
needs the parameters of a complex interval: xx is a real interval
from which a surrounding complex interval is constructed. r (>1)
controls its size. In the current version, r is set to 4.0 in
cinte.p
- The second function fb1
needs a real interval xx as input. Then a b_interval is
constructed containing additionally a bit that gives information
about boundedness.
- The third function f1
is a simple interval function that is needed in the final phase
of the algorithm, where quadrature formulas are applied.
begin
....
init_cinte; { preparation for integration - reads all coefficients }
{ of the used Gaussian quadrature formulas }
....
( it is faster to read them only once if the integration routine is
called more than once )
....
integral:= integrate(f1,fb1,fa1, 10.0, 40.0, 1.0e-12);
....
calculates the integral
with an error (see the description of eps above) <1.0e-12
Compilation:
(make if necessary)
mxsc test1 . The executable is the file test1.
Use of cinte_comf.p:
This program is simpler to use but slower than cinte.p
(In the average of the tests, cinte_comf.p was about 20%
slower than cinte.p)
There are two main differences:
- There is only one function with argument and result type
int_type
- You do not have construct the arguments within the function.
Our example therefore reads
program test (input, output);
use cinte_comf;
....
function f1(x : inte_type): inte_type;
begin
f1 := sin(sqr(x))*exp(x);
end;
....
begin
cinte_init;
....
integral:= integrate(f1, 10.0, 40.0, 1.0e-12);
....
Compilation:
(make if necessary)
mxsc test3 . The executable is the file test3.
Data types in cinte:
In order to calculate a bound for the integrand we use the module
bi_ari. In this module, the type
global type b_interval = global record
int : interval;
b : boolean;
end;
is defined and an arithmetic is programmed. It is essentially interval
arithmetic but the boolean value controls the boundedness.
Suppose an operation has to be applied to a b_interval.
If b=1
(the initial value), we check if the next operation may yield
a floating point exception. If there will be no floating point exception,
the operation is applied to the interval
and b=1 is kept. Otherwise, the interval is set to [0,0] and the boolean to 0.
If b=0, we return the input value.
In order to calculate a bound for the function in a complex neighborhood of
an interval, we use the module
ade_ari. In this module, the type
global type analyticd = global record
rint, dint : interval;
cint : icom;
a, d : boolean;
end;
is defined. rint and dint are bounds for the values
and first derivatives respectively. cint is a complex interval.
The boolean variables a,d control analyticity of a function/operation
on cint and differentiability on rint. Derivatives are
calculated with Automatic Differentiation.
The type analyticd contains a component of type icom
This represents a complex (rectangular) interval that is symmetric
with respect to the real axis:
global type icom = global record
re : interval;
im : real;
end;
It is particular designed for real problems, i.e., with problems that
deal with functions that are analytic as well as real-valued for
real arguments. The corresponding arithmetic is faster than
arithmetic for the type cinterval
Data type in cinte_comf:
The type
global type inte_type analyticd = global record
rint, dint : interval;
cint : icom;
a, d : boolean;
end;
(which is essentially the same as analyticd)
serves as a container for interval (only rint is used),
for b_interval (only rint and a are used
analogously to the components int and b in
b_interval) and for analyticd.
Switch between the three types represented by inte_type
is done in cinte_comf.p
by the global variable type_switch
(1=interval-operations, 2=b_interval-operations, 3=analyticd-operations)
Available functions:
At the moment, the available functions / operations for argument(s)
of the type analyticd (or inte_type) are
- +,-,*, / : combinations of analyticd
(or inte_type)
with analyticd
(or inte_type resp.),
interval, real, integer
- write
- construct (in ade_ari.p) :
Arguments: (x:interval;r:real) or
(x:interval;r1,r2:real); constructs analyticd that consists of a real interval
...rint=x, a derivative ...dint=1 and a complex interval ...cint
with real part
mid(x)+(x-mid(x))*r1 and imaginary part (x-mid(x))*r2. If only r is
given, we set r1 to (r+1/r)/2 and r2 to (r-1/r)/2.
Differentiability ...d and analyticity ...a are set to TRUE.
- iconstr (in ade_comf.p) :
same as construct, but yields an
inte_type and additionally
to the previous
description, we can use also only the argument (x:interval).
It initializes ...rint=x and ...a=TRUE (representing an interval x
and the knowledge of its boundedness TRUE).
- - (sign change), sqr, exp,
sin, cos
- pow(x,n) : yields the nth power (n: integer) of x
- pospart : the positive part extended to an
analytic function wherever possible.
- abs : the absolute value extended to an
analytic function wherever possible.
- sqrtabs :
square root of the absolute value analytically extended to
an
analytic function wherever possible
In arith.p, these functions, except for pow,
are also implemented for the type
icom. Furthermore, tan and ln
are provided.