# calculus - Man Page

Integration and ordinary differential equations

## Synopsis

`package require `

**Tcl 8.4**

`package require `

**math::calculus 0.8.2**

**::math::calculus::integral** *begin end nosteps func*

**::math::calculus::integralExpr** *begin end nosteps expression*

**::math::calculus::integral2D** *xinterval yinterval func*

**::math::calculus::integral2D_accurate** *xinterval yinterval func*

**::math::calculus::integral3D** *xinterval yinterval zinterval func*

**::math::calculus::integral3D_accurate** *xinterval yinterval zinterval func*

**::math::calculus::qk15** *xstart xend func nosteps*

**::math::calculus::qk15_detailed** *xstart xend func nosteps*

**::math::calculus::eulerStep** *t tstep xvec func*

**::math::calculus::heunStep** *t tstep xvec func*

**::math::calculus::rungeKuttaStep** *t tstep xvec func*

**::math::calculus::boundaryValueSecondOrder** *coeff_func force_func leftbnd rightbnd nostep*

**::math::calculus::solveTriDiagonal** *acoeff bcoeff ccoeff dvalue*

**::math::calculus::newtonRaphson** *func deriv initval*

**::math::calculus::newtonRaphsonParameters** *maxiter tolerance*

**::math::calculus::regula_falsi** *f xb xe eps*

## Description

This package implements several simple mathematical algorithms:

- The integration of a function over an interval
- The numerical integration of a system of ordinary differential equations.
- Estimating the root(s) of an equation of one variable.

The package is fully implemented in Tcl. No particular attention has been paid to the accuracy of the calculations. Instead, well-known algorithms have been used in a straightforward manner.

This document describes the procedures and explains their usage.

## Procedures

This package defines the following public procedures:

**::math::calculus::integral***begin end nosteps func*Determine the integral of the given function using the Simpson rule. The interval for the integration is [

*begin*,*end*]. The remaining arguments are:*nosteps*Number of steps in which the interval is divided.

*func*Function to be integrated. It should take one single argument.

**::math::calculus::integralExpr***begin end nosteps expression*Similar to the previous proc, this one determines the integral of the given

*expression*using the Simpson rule. The interval for the integration is [*begin*,*end*]. The remaining arguments are:*nosteps*Number of steps in which the interval is divided.

*expression*Expression to be integrated. It should use the variable "x" as the only variable (the "integrate")

**::math::calculus::integral2D***xinterval yinterval func***::math::calculus::integral2D_accurate***xinterval yinterval func*The commands

**integral2D**and**integral2D_accurate**calculate the integral of a function of two variables over the rectangle given by the first two arguments, each a list of three items, the start and stop interval for the variable and the number of steps.The command

**integral2D**evaluates the function at the centre of each rectangle, whereas the command**integral2D_accurate**uses a four-point quadrature formula. This results in an exact integration of polynomials of third degree or less.The function must take two arguments and return the function value.

**::math::calculus::integral3D***xinterval yinterval zinterval func***::math::calculus::integral3D_accurate***xinterval yinterval zinterval func*The commands

**integral3D**and**integral3D_accurate**are the three-dimensional equivalent of**integral2D**and**integral3D_accurate**. The function*func*takes three arguments and is integrated over the block in 3D space given by three intervals.**::math::calculus::qk15***xstart xend func nosteps*Determine the integral of the given function using the Gauss-Kronrod 15 points quadrature rule. The returned value is the estimate of the integral over the interval [

*xstart*,*xend*]. The remaining arguments are:*func*Function to be integrated. It should take one single argument.

- ?nosteps?
Number of steps in which the interval is divided. Defaults to 1.

**::math::calculus::qk15_detailed***xstart xend func nosteps*Determine the integral of the given function using the Gauss-Kronrod 15 points quadrature rule. The interval for the integration is [

*xstart*,*xend*]. The procedure returns a list of four values:- The estimate of the integral over the specified interval (I).
- An estimate of the absolute error in I.
- The estimate of the integral of the absolute value of the function over the interval.
- The estimate of the integral of the absolute value of the function minus its mean over the interval.

The remaining arguments are:

*func*Function to be integrated. It should take one single argument.

- ?nosteps?
Number of steps in which the interval is divided. Defaults to 1.

**::math::calculus::eulerStep***t tstep xvec func*Set a single step in the numerical integration of a system of differential equations. The method used is Euler's.

*t*Value of the independent variable (typically time) at the beginning of the step.

*tstep*Step size for the independent variable.

*xvec*List (vector) of dependent values

*func*Function of t and the dependent values, returning a list of the derivatives of the dependent values. (The lengths of xvec and the return value of "func" must match).

**::math::calculus::heunStep***t tstep xvec func*Set a single step in the numerical integration of a system of differential equations. The method used is Heun's.

*t*Value of the independent variable (typically time) at the beginning of the step.

*tstep*Step size for the independent variable.

*xvec*List (vector) of dependent values

*func*Function of t and the dependent values, returning a list of the derivatives of the dependent values. (The lengths of xvec and the return value of "func" must match).

**::math::calculus::rungeKuttaStep***t tstep xvec func*Set a single step in the numerical integration of a system of differential equations. The method used is Runge-Kutta 4th order.

*t*Value of the independent variable (typically time) at the beginning of the step.

*tstep*Step size for the independent variable.

*xvec*List (vector) of dependent values

*func*Function of t and the dependent values, returning a list of the derivatives of the dependent values. (The lengths of xvec and the return value of "func" must match).

**::math::calculus::boundaryValueSecondOrder***coeff_func force_func leftbnd rightbnd nostep*Solve a second order linear differential equation with boundary values at two sides. The equation has to be of the form (the "conservative" form):

d dy d -- A(x)-- + -- B(x)y + C(x)y = D(x) dx dx dx

Ordinarily, such an equation would be written as:

d2y dy a(x)--- + b(x)-- + c(x) y = D(x) dx2 dx

The first form is easier to discretise (by integrating over a finite volume) than the second form. The relation between the two forms is fairly straightforward:

A(x) = a(x) B(x) = b(x) - a'(x) C(x) = c(x) - B'(x) = c(x) - b'(x) + a''(x)

Because of the differentiation, however, it is much easier to ask the user to provide the functions A, B and C directly.

*coeff_func*Procedure returning the three coefficients (A, B, C) of the equation, taking as its one argument the x-coordinate.

*force_func*Procedure returning the right-hand side (D) as a function of the x-coordinate.

*leftbnd*A list of two values: the x-coordinate of the left boundary and the value at that boundary.

*rightbnd*A list of two values: the x-coordinate of the right boundary and the value at that boundary.

*nostep*Number of steps by which to discretise the interval. The procedure returns a list of x-coordinates and the approximated values of the solution.

**::math::calculus::solveTriDiagonal***acoeff bcoeff ccoeff dvalue*Solve a system of linear equations Ax = b with A a tridiagonal matrix. Returns the solution as a list.

*acoeff*List of values on the lower diagonal

*bcoeff*List of values on the main diagonal

*ccoeff*List of values on the upper diagonal

*dvalue*List of values on the righthand-side

**::math::calculus::newtonRaphson***func deriv initval*Determine the root of an equation given by

func(x) = 0

using the method of Newton-Raphson. The procedure takes the following arguments:

*func*Procedure that returns the value the function at x

*deriv*Procedure that returns the derivative of the function at x

*initval*Initial value for x

**::math::calculus::newtonRaphsonParameters***maxiter tolerance*Set the numerical parameters for the Newton-Raphson method:

*maxiter*Maximum number of iteration steps (defaults to 20)

*tolerance*Relative precision (defaults to 0.001)

**::math::calculus::regula_falsi***f xb xe eps*Return an estimate of the zero or one of the zeros of the function contained in the interval [xb,xe]. The error in this estimate is of the order of eps*abs(xe-xb), the actual error may be slightly larger.

The method used is the so-called

*regula falsi*or*false position*method. It is a straightforward implementation. The method is robust, but requires that the interval brackets a zero or at least an uneven number of zeros, so that the value of the function at the start has a different sign than the value at the end.In contrast to Newton-Raphson there is no need for the computation of the function's derivative.

- command
*f* Name of the command that evaluates the function for which the zero is to be returned

- float
*xb* Start of the interval in which the zero is supposed to lie

- float
*xe* End of the interval

- float
*eps* Relative allowed error (defaults to 1.0e-4)

- command

*Notes:*

Several of the above procedures take the *names* of procedures as arguments. To avoid problems with the *visibility* of these procedures, the fully-qualified name of these procedures is determined inside the calculus routines. For the user this has only one consequence: the named procedure must be visible in the calling procedure. For instance:

namespace eval ::mySpace { namespace export calcfunc proc calcfunc { x } { return $x } } # # Use a fully-qualified name # namespace eval ::myCalc { proc detIntegral { begin end } { return [integral $begin $end 100 ::mySpace::calcfunc] } } # # Import the name # namespace eval ::myCalc { namespace import ::mySpace::calcfunc proc detIntegral { begin end } { return [integral $begin $end 100 calcfunc] } }

Enhancements for the second-order boundary value problem:

- Other types of boundary conditions (zero gradient, zero flux)
- Other schematisation of the first-order term (now central differences are used, but upstream differences might be useful too).

## Examples

Let us take a few simple examples:

Integrate x over the interval [0,100] (20 steps):

proc linear_func { x } { return $x } puts "Integral: [::math::calculus::integral 0 100 20 linear_func]"

For simple functions, the alternative could be:

puts "Integral: [::math::calculus::integralExpr 0 100 20 {$x}]"

Do not forget the braces!

The differential equation for a dampened oscillator:

x'' + rx' + wx = 0

can be split into a system of first-order equations:

x' = y y' = -ry - wx

Then this system can be solved with code like this:

proc dampened_oscillator { t xvec } { set x [lindex $xvec 0] set x1 [lindex $xvec 1] return [list $x1 [expr {-$x1-$x}]] } set xvec { 1.0 0.0 } set t 0.0 set tstep 0.1 for { set i 0 } { $i < 20 } { incr i } { set result [::math::calculus::eulerStep $t $tstep $xvec dampened_oscillator] puts "Result ($t): $result" set t [expr {$t+$tstep}] set xvec $result }

Suppose we have the boundary value problem:

Dy'' + ky = 0 x = 0: y = 1 x = L: y = 0

This boundary value problem could originate from the diffusion of a decaying substance.

It can be solved with the following fragment:

proc coeffs { x } { return [list $::Diff 0.0 $::decay] } proc force { x } { return 0.0 } set Diff 1.0e-2 set decay 0.0001 set length 100.0 set y [::math::calculus::boundaryValueSecondOrder \ coeffs force {0.0 1.0} [list $length 0.0] 100]

## Bugs, Ideas, Feedback

This document, and the package it describes, will undoubtedly contain bugs and other problems. Please report such in the category *math :: calculus* of the *Tcllib Trackers* [http://core.tcl.tk/tcllib/reportlist]. Please also report any ideas for enhancements you may have for either package and/or documentation.

When proposing code changes, please provide *unified diffs*, i.e the output of **diff -u**.

Note further that *attachments* are strongly preferred over inlined patches. Attachments can be made by going to the **Edit** form of the ticket immediately after its creation, and then using the left-most button in the secondary navigation bar.

## See Also

romberg

## Keywords

calculus, differential equations, integration, math, roots

## Category

Mathematics

## Copyright

Copyright (c) 2002,2003,2004 Arjen Markus