math::calculus(n)

📦

version 2.0 · 2 versions

2.0 · fedora 43, fedora 44, nixos 26.05, nixpkgs unstable (latest)
1.21 · alpine 3.24, epel 9, opensuse-leap 16.0

, in fedora 43, fedora 44, nixos 26.05, nixpkgs unstable

✨ math::calculus — the Tcl Math Library ✨

okay bestie, grab a snack, we're doing math today and honestly? it's giving iconic 💅

🎀 NAME 🎀

math::calculus — integration and ordinary differential equations. That's the whole vibe. 🧮✨

🌸 SYNOPSIS 🌸

package require Tcl 8.5 9

package require math::calculus 1.1

::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

::math::calculus::root_bisection f xb xe eps

::math::calculus::root_secant f xb xe eps

::math::calculus::root_brent f xb xe eps

::math::calculus::root_chandrupatla f xb xe eps

💖 DESCRIPTION 💖

This package is a lil collection of simple mathematical algorithms, and here's the lineup:

Integrating 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

and yes, it's fully implemented in Tcl — no C, no compiling, just vibes and pure Tcl, bestie. 💅

⚠️ Real talk though: no particular attention has been paid to the accuracy of the calculations. Instead, well-known algorithms have just been used in a straightforward manner. So like… she works, but she is NOT out here chasing the last decimal. Manage your expectations, fr fr. 😭

This document describes the procedures and explains their usage. Slay.

👑 PROCEDURES 👑

This package defines the following public procedures:

::math::calculus::integral begin end nosteps func — Figure out 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 — the function to be integrated. It should take one single argument. One. Just the one, okay? 🙈
::math::calculus::integralExpr begin end nosteps expression — Same energy as the previous proc, except this one works out 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 — the 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

So integral2D and integral2D_accurate both calculate the integral of a function of two variables over the rectangle given by the first two arguments — and each of those arguments is a list of three items: the start and stop interval for the variable, plus the number of steps. 📋

integral2D evaluates the function at the centre of each rectangle, whereas integral2D_accurate uses a four-point quadrature formula. Which means… integral2D_accurate gives you an exact integration of polynomials of third degree or less. No cap. 🔥

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

These two are the three-dimensional equivalent of integral2D and integral3D_accurate. The function func takes three arguments and gets integrated over the block in 3D space given by three intervals. It's giving… volume. 🦋

::math::calculus::qk15 xstart xend func nosteps — Determine the integral of the given function using the Gauss-Kronrod 15 points quadrature rule</cornflower>. The returned value is the estimate of the integral over the interval [xstart, xend]. The remaining arguments are:
- func — the function to be integrated. It should take one single argument.
- ?nosteps? — number of steps in which the interval is divided. Defaults to 1. (Note the lil question marks — she's optional, bestie. 🫶)
::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 hands back a list of four values:
1. The estimate of the integral over the specified interval (I).
2. An estimate of the absolute error in I.
3. The estimate of the integral of the absolute value of the function over the interval.
4. The estimate of the integral of the absolute value of the function minus its mean over the interval.

The remaining arguments are:

func — the function to be integrated. It should take one single argument.
?nosteps? — number of steps in which the interval is divided. Defaults to 1. (Optional again, iykyk 💅)

::math::calculus::eulerStep t tstep xvec func — Take 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 — like, they HAVE to match, okay? 💀)
::math::calculus::heunStep t tstep xvec func — Take 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 — Take a single step in the numerical integration of a system of differential equations. The method used is Runge-Kutta 4th order. She's the fancy one and honestly? deserved. 👑
- 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 all that differentiation, though, it's much easier to just ask the user to hand over the functions A, B and C directly. So that's what we do. 💅

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 where A is 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) — and heads up bestie, the actual error may be slightly larger. 🙈

okay that's the crew so far, but we're not done yet — more parts are coming, so stay tuned 💖

okay so bestie, the method used here 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. periodt. ✨

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) 💅

`::math::calculus::root_bisection 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 bisection. For properties: see the regula falsi procedure. 🦋

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-7) 💅

`::math::calculus::root_secant f xb xe eps`

The method used is the so-called secant method. For properties: see the regula falsi procedure. Note that this method is not guaranteed to produce a zero, but it is fast. 🔥

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-7) 💅

`::math::calculus::root_brent f xb xe eps`

The method used is the so-called Brent method. For properties: see the regula falsi procedure. 🌈

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-7) 💅

`::math::calculus::root_chandrupatla f xb xe eps`

The method used is the so-called Chandrupatla method. It is a variant of the Brent method and is reputed to be faster. See https://www.embeddedrelated.com/showarticle/855.php for a discussion. 👑

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-7) 💅

📝 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). 🌀
The various root finding methods differ in robustness, ease of use and convergence rates. While the Newton-Raphson method is quite fast (quadratic convergence), it is not guaranteed to produce an answer. The methods that bracket the root, require two starting points, but except for the secant method are guaranteed to deliver a good estimate. 💯

✨ 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

that's the whole page, bestie — from tcllib 1.1, the math::calculus(n) page. no date given, just vibes. 💖✨

📦 where this page comes from

fedora 43: tcllib 2.0-3.fc43 · TCL
fedora 44: tcllib 2.0-4.fc44 · TCL
nixos 26.05: tcllib 2.0
nixpkgs unstable: tcllib 2.0
💖 girlified by deepseek-ai/DeepSeek-V4.1-Flash on 2026-10-18: a language model's rewrite of the original, under its terms