=============================================================================

	icalc - a complex-number expression language

	(C) 1991, 1992 Martin W. Scott. All Rights Reserved.

=============================================================================

	Integration Notes

=============================================================================

There are a few icalc scripts that deal with the problem of real
integration. These scripts provide 4 basic routines which vary in their
sophistication, and also in their speed of calculation. They are:

	gl.ic		- contains definition of function gl, which
			  calculates integrals using a 10-point
			  Gauss-Legendre quadrature.

	trapzd.ic	- contains definition of function trapzd, used
			  by qtrap, qsimp and qromb routines.

	polint.ic	- contains definition of function polint, used
			  by qromb, but useful in its own right.

	qtrap.ic	- contains definition of function qtrap, which
			  performs numerical integration using the
			  Trapezium Rule.

	qsimp.ic	- contains definition of function qtrap, which
			  performs numerical integration using the
			  Simpson's Rule.

	qromb.ic	- contains definition of function qtrap, which
			  performs Romberg numerical integration 

All the routines assume that the array base is 1, i.e. given an array a,
its first element is a[1]. If you have changed the array base (with a
call to arraybase()) you should restore its value to 1 before using
these routines.



How to use the integration routines
-----------------------------------
Start icalc with the desired routine, from the Workbench or a Shell.
Here, xxx denotes gl, qromb, qtrap or qsimp.

	From Shell:	type "icalc xxx.ic -" at shell prompt

	From WB:	click once on icalc icon, hold down shift
			and double-click on xxx.ic icon.

All the routines gl, qtrap, qsimp and qromb take 3 arguments - an
expression in x, and the interval endpoints a and b.

The expression in x can be any expression at all; here are some
examples:

	> gl(x*x,0,1)
		0.33333333328

	> gl(sqr(x),0,1)
		0.33333333328

	> gl(sin(ln(x)),1,2)
		0.36972237492

	> qtrap(sin(ln(x)),1,2)
		0.36972159245

	> qsimp(log(cos(x)),0,PI/4)
		-0.03752900436

	> 2*sqr(qromb(exp(-sqr(x)/2), 0, 100))
		3.14159314113

and so on. The question remains, when to use which routine?

The qtrap routine may be used for integrating functions which are not
very smooth; since most expression-based functions are, you should
hardly ever need to use it. It is the slowest of the methods, and can
take many seconds to produce a result.

The qsimp routine is usually more efficient than qtrap when the
integrand has a continuous 3rd derivative.

The qromb routine is best for sufficiently smooth functions (e.g.
analytic functions) and is usually fastest at producing a result.

The gl routine is useful for quick results, although usually qromb is
quick enough (about a second or two) for most integrations. gl is useful
when used in a routine for multi-dimensional integration.

For more information on all routines, consult "Numerical Recipes in C",
by Press et al. That is the source of these algorithms.



Improving Accuracy
------------------
You can improve the accuracy of results by setting the global variables
controlling whichever routine you are using (this does not apply to gl).
The variables, with their default values, are as follows:

	Routine		Variable	Value
	-------		--------	-----
	qtrap		QTRAP_EPS	1e-5
	qsimp		QSIMP_EPS	1e-6
	qromb		QROMB_EPS	1e-6

The xxx_EPS variables denote the fractional accuracy desired, i.e. if
the (error estimate ds) < xxx_EPS*(integral estimate s) then the routine
returns the integral estimate, s.

In some cases, desired accuracy cannot be achieved within the maximum
allowed iterations, causing a user error. In such cases, try another
routine.

Finally, a warning: because of the type of termination condition
(fractional accuracy) integrations with a result of zero cause problems.
If your routine is taking an excessively long time (more than 10 seconds
for qromb, say) it may be because the result is zero. Try splitting the
integration limits up, and doing two integrations. If this fails, graph
the function to see what is going on. The routines may also have
problems with functions containing singularities.

Examples (all using best routine, qromb):

	> qromb(sin(x),0,PI)
		1.99999999459
	> QROMB_EPS=1e-8
		0.00000001
	> qromb(sin(x),0,PI)
		2

	> qromb(sin(x),0,2*PI)

No result appears for some time. Press control-c to abort.

	<CONTROL>-C
	*** break
	> qromb(sin(x),0,PI) + qromb(sin(x),PI,2*PI)
		-4.21884749358e-15

The last integration was very close to zero, which is why the
qromb(sin(x),0,2*PI) call was taking so long.

That's about all I want to say about these routines. The scripts
containing them are partially-commented, but you'd best get a hold of
Numerical Recipes if you need more information.