'\" '\" Generated from file 'calculus\&.man' by tcllib/doctools with format 'nroff' '\" Copyright (c) 2002,2003,2004 Arjen Markus '\" .TH "math::calculus" 3tcl 0\&.8\&.1 tcllib "Tcl Math Library" .\" The -*- nroff -*- definitions below are for supplemental macros used .\" in Tcl/Tk manual entries. .\" .\" .AP type name in/out ?indent? .\" Start paragraph describing an argument to a library procedure. .\" type is type of argument (int, etc.), in/out is either "in", "out", .\" or "in/out" to describe whether procedure reads or modifies arg, .\" and indent is equivalent to second arg of .IP (shouldn't ever be .\" needed; use .AS below instead) .\" .\" .AS ?type? ?name? .\" Give maximum sizes of arguments for setting tab stops. Type and .\" name are examples of largest possible arguments that will be passed .\" to .AP later. If args are omitted, default tab stops are used. .\" .\" .BS .\" Start box enclosure. From here until next .BE, everything will be .\" enclosed in one large box. .\" .\" .BE .\" End of box enclosure. .\" .\" .CS .\" Begin code excerpt. .\" .\" .CE .\" End code excerpt. .\" .\" .VS ?version? ?br? .\" Begin vertical sidebar, for use in marking newly-changed parts .\" of man pages. The first argument is ignored and used for recording .\" the version when the .VS was added, so that the sidebars can be .\" found and removed when they reach a certain age. If another argument .\" is present, then a line break is forced before starting the sidebar. .\" .\" .VE .\" End of vertical sidebar. .\" .\" .DS .\" Begin an indented unfilled display. .\" .\" .DE .\" End of indented unfilled display. .\" .\" .SO ?manpage? .\" Start of list of standard options for a Tk widget. The manpage .\" argument defines where to look up the standard options; if .\" omitted, defaults to "options". The options follow on successive .\" lines, in three columns separated by tabs. .\" .\" .SE .\" End of list of standard options for a Tk widget. .\" .\" .OP cmdName dbName dbClass .\" Start of description of a specific option. cmdName gives the .\" option's name as specified in the class command, dbName gives .\" the option's name in the option database, and dbClass gives .\" the option's class in the option database. .\" .\" .UL arg1 arg2 .\" Print arg1 underlined, then print arg2 normally. .\" .\" .QW arg1 ?arg2? .\" Print arg1 in quotes, then arg2 normally (for trailing punctuation). .\" .\" .PQ arg1 ?arg2? .\" Print an open parenthesis, arg1 in quotes, then arg2 normally .\" (for trailing punctuation) and then a closing parenthesis. .\" .\" # Set up traps and other miscellaneous stuff for Tcl/Tk man pages. .if t .wh -1.3i ^B .nr ^l \n(.l .ad b .\" # Start an argument description .de AP .ie !"\\$4"" .TP \\$4 .el \{\ . ie !"\\$2"" .TP \\n()Cu . el .TP 15 .\} .ta \\n()Au \\n()Bu .ie !"\\$3"" \{\ \&\\$1 \\fI\\$2\\fP (\\$3) .\".b .\} .el \{\ .br .ie !"\\$2"" \{\ \&\\$1 \\fI\\$2\\fP .\} .el \{\ \&\\fI\\$1\\fP .\} .\} .. .\" # define tabbing values for .AP .de AS .nr )A 10n .if !"\\$1"" .nr )A \\w'\\$1'u+3n .nr )B \\n()Au+15n .\" .if !"\\$2"" .nr )B \\w'\\$2'u+\\n()Au+3n .nr )C \\n()Bu+\\w'(in/out)'u+2n .. .AS Tcl_Interp Tcl_CreateInterp in/out .\" # BS - start boxed text .\" # ^y = starting y location .\" # ^b = 1 .de BS .br .mk ^y .nr ^b 1u .if n .nf .if n .ti 0 .if n \l'\\n(.lu\(ul' .if n .fi .. .\" # BE - end boxed text (draw box now) .de BE .nf .ti 0 .mk ^t .ie n \l'\\n(^lu\(ul' .el \{\ .\" Draw four-sided box normally, but don't draw top of .\" box if the box started on an earlier page. .ie !\\n(^b-1 \{\ \h'-1.5n'\L'|\\n(^yu-1v'\l'\\n(^lu+3n\(ul'\L'\\n(^tu+1v-\\n(^yu'\l'|0u-1.5n\(ul' .\} .el \}\ \h'-1.5n'\L'|\\n(^yu-1v'\h'\\n(^lu+3n'\L'\\n(^tu+1v-\\n(^yu'\l'|0u-1.5n\(ul' .\} .\} .fi .br .nr ^b 0 .. .\" # VS - start vertical sidebar .\" # ^Y = starting y location .\" # ^v = 1 (for troff; for nroff this doesn't matter) .de VS .if !"\\$2"" .br .mk ^Y .ie n 'mc \s12\(br\s0 .el .nr ^v 1u .. .\" # VE - end of vertical sidebar .de VE .ie n 'mc .el \{\ .ev 2 .nf .ti 0 .mk ^t \h'|\\n(^lu+3n'\L'|\\n(^Yu-1v\(bv'\v'\\n(^tu+1v-\\n(^Yu'\h'-|\\n(^lu+3n' .sp -1 .fi .ev .\} .nr ^v 0 .. .\" # Special macro to handle page bottom: finish off current .\" # box/sidebar if in box/sidebar mode, then invoked standard .\" # page bottom macro. .de ^B .ev 2 'ti 0 'nf .mk ^t .if \\n(^b \{\ .\" Draw three-sided box if this is the box's first page, .\" draw two sides but no top otherwise. .ie !\\n(^b-1 \h'-1.5n'\L'|\\n(^yu-1v'\l'\\n(^lu+3n\(ul'\L'\\n(^tu+1v-\\n(^yu'\h'|0u'\c .el \h'-1.5n'\L'|\\n(^yu-1v'\h'\\n(^lu+3n'\L'\\n(^tu+1v-\\n(^yu'\h'|0u'\c .\} .if \\n(^v \{\ .nr ^x \\n(^tu+1v-\\n(^Yu \kx\h'-\\nxu'\h'|\\n(^lu+3n'\ky\L'-\\n(^xu'\v'\\n(^xu'\h'|0u'\c .\} .bp 'fi .ev .if \\n(^b \{\ .mk ^y .nr ^b 2 .\} .if \\n(^v \{\ .mk ^Y .\} .. .\" # DS - begin display .de DS .RS .nf .sp .. .\" # DE - end display .de DE .fi .RE .sp .. .\" # SO - start of list of standard options .de SO 'ie '\\$1'' .ds So \\fBoptions\\fR 'el .ds So \\fB\\$1\\fR .SH "STANDARD OPTIONS" .LP .nf .ta 5.5c 11c .ft B .. .\" # SE - end of list of standard options .de SE .fi .ft R .LP See the \\*(So manual entry for details on the standard options. .. .\" # OP - start of full description for a single option .de OP .LP .nf .ta 4c Command-Line Name: \\fB\\$1\\fR Database Name: \\fB\\$2\\fR Database Class: \\fB\\$3\\fR .fi .IP .. .\" # CS - begin code excerpt .de CS .RS .nf .ta .25i .5i .75i 1i .. .\" # CE - end code excerpt .de CE .fi .RE .. .\" # UL - underline word .de UL \\$1\l'|0\(ul'\\$2 .. .\" # QW - apply quotation marks to word .de QW .ie '\\*(lq'"' ``\\$1''\\$2 .\"" fix emacs highlighting .el \\*(lq\\$1\\*(rq\\$2 .. .\" # PQ - apply parens and quotation marks to word .de PQ .ie '\\*(lq'"' (``\\$1''\\$2)\\$3 .\"" fix emacs highlighting .el (\\*(lq\\$1\\*(rq\\$2)\\$3 .. .\" # QR - quoted range .de QR .ie '\\*(lq'"' ``\\$1''\\-``\\$2''\\$3 .\"" fix emacs highlighting .el \\*(lq\\$1\\*(rq\\-\\*(lq\\$2\\*(rq\\$3 .. .\" # MT - "empty" string .de MT .QW "" .. .BS .SH NAME math::calculus \- Integration and ordinary differential equations .SH SYNOPSIS package require \fBTcl 8\&.4\fR .sp package require \fBmath::calculus 0\&.8\&.1\fR .sp \fB::math::calculus::integral\fR \fIbegin\fR \fIend\fR \fInosteps\fR \fIfunc\fR .sp \fB::math::calculus::integralExpr\fR \fIbegin\fR \fIend\fR \fInosteps\fR \fIexpression\fR .sp \fB::math::calculus::integral2D\fR \fIxinterval\fR \fIyinterval\fR \fIfunc\fR .sp \fB::math::calculus::integral2D_accurate\fR \fIxinterval\fR \fIyinterval\fR \fIfunc\fR .sp \fB::math::calculus::integral3D\fR \fIxinterval\fR \fIyinterval\fR \fIzinterval\fR \fIfunc\fR .sp \fB::math::calculus::integral3D_accurate\fR \fIxinterval\fR \fIyinterval\fR \fIzinterval\fR \fIfunc\fR .sp \fB::math::calculus::qk15\fR \fIxstart\fR \fIxend\fR \fIfunc\fR \fInosteps\fR .sp \fB::math::calculus::qk15_detailed\fR \fIxstart\fR \fIxend\fR \fIfunc\fR \fInosteps\fR .sp \fB::math::calculus::eulerStep\fR \fIt\fR \fItstep\fR \fIxvec\fR \fIfunc\fR .sp \fB::math::calculus::heunStep\fR \fIt\fR \fItstep\fR \fIxvec\fR \fIfunc\fR .sp \fB::math::calculus::rungeKuttaStep\fR \fIt\fR \fItstep\fR \fIxvec\fR \fIfunc\fR .sp \fB::math::calculus::boundaryValueSecondOrder\fR \fIcoeff_func\fR \fIforce_func\fR \fIleftbnd\fR \fIrightbnd\fR \fInostep\fR .sp \fB::math::calculus::solveTriDiagonal\fR \fIacoeff\fR \fIbcoeff\fR \fIccoeff\fR \fIdvalue\fR .sp \fB::math::calculus::newtonRaphson\fR \fIfunc\fR \fIderiv\fR \fIinitval\fR .sp \fB::math::calculus::newtonRaphsonParameters\fR \fImaxiter\fR \fItolerance\fR .sp \fB::math::calculus::regula_falsi\fR \fIf\fR \fIxb\fR \fIxe\fR \fIeps\fR .sp .BE .SH DESCRIPTION .PP This package implements several simple mathematical algorithms: .IP \(bu The integration of a function over an interval .IP \(bu The numerical integration of a system of ordinary differential equations\&. .IP \(bu Estimating the root(s) of an equation of one variable\&. .PP .PP 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\&. .PP This document describes the procedures and explains their usage\&. .SH PROCEDURES This package defines the following public procedures: .TP \fB::math::calculus::integral\fR \fIbegin\fR \fIend\fR \fInosteps\fR \fIfunc\fR Determine the integral of the given function using the Simpson rule\&. The interval for the integration is [\fIbegin\fR, \fIend\fR]\&. The remaining arguments are: .RS .TP \fInosteps\fR Number of steps in which the interval is divided\&. .TP \fIfunc\fR Function to be integrated\&. It should take one single argument\&. .RE .sp .TP \fB::math::calculus::integralExpr\fR \fIbegin\fR \fIend\fR \fInosteps\fR \fIexpression\fR Similar to the previous proc, this one determines the integral of the given \fIexpression\fR using the Simpson rule\&. The interval for the integration is [\fIbegin\fR, \fIend\fR]\&. The remaining arguments are: .RS .TP \fInosteps\fR Number of steps in which the interval is divided\&. .TP \fIexpression\fR Expression to be integrated\&. It should use the variable "x" as the only variable (the "integrate") .RE .sp .TP \fB::math::calculus::integral2D\fR \fIxinterval\fR \fIyinterval\fR \fIfunc\fR .TP \fB::math::calculus::integral2D_accurate\fR \fIxinterval\fR \fIyinterval\fR \fIfunc\fR The commands \fBintegral2D\fR and \fBintegral2D_accurate\fR 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\&. .sp The command \fBintegral2D\fR evaluates the function at the centre of each rectangle, whereas the command \fBintegral2D_accurate\fR uses a four-point quadrature formula\&. This results in an exact integration of polynomials of third degree or less\&. .sp The function must take two arguments and return the function value\&. .TP \fB::math::calculus::integral3D\fR \fIxinterval\fR \fIyinterval\fR \fIzinterval\fR \fIfunc\fR .TP \fB::math::calculus::integral3D_accurate\fR \fIxinterval\fR \fIyinterval\fR \fIzinterval\fR \fIfunc\fR The commands \fBintegral3D\fR and \fBintegral3D_accurate\fR are the three-dimensional equivalent of \fBintegral2D\fR and \fBintegral3D_accurate\fR\&. The function \fIfunc\fR takes three arguments and is integrated over the block in 3D space given by three intervals\&. .TP \fB::math::calculus::qk15\fR \fIxstart\fR \fIxend\fR \fIfunc\fR \fInosteps\fR 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 [\fIxstart\fR, \fIxend\fR]\&. The remaining arguments are: .RS .TP \fIfunc\fR Function to be integrated\&. It should take one single argument\&. .TP ?nosteps? Number of steps in which the interval is divided\&. Defaults to 1\&. .RE .sp .TP \fB::math::calculus::qk15_detailed\fR \fIxstart\fR \fIxend\fR \fIfunc\fR \fInosteps\fR Determine the integral of the given function using the Gauss-Kronrod 15 points quadrature rule\&. The interval for the integration is [\fIxstart\fR, \fIxend\fR]\&. The procedure returns a list of four values: .RS .IP \(bu The estimate of the integral over the specified interval (I)\&. .IP \(bu An estimate of the absolute error in I\&. .IP \(bu The estimate of the integral of the absolute value of the function over the interval\&. .IP \(bu The estimate of the integral of the absolute value of the function minus its mean over the interval\&. .RE .IP The remaining arguments are: .RS .TP \fIfunc\fR Function to be integrated\&. It should take one single argument\&. .TP ?nosteps? Number of steps in which the interval is divided\&. Defaults to 1\&. .RE .sp .TP \fB::math::calculus::eulerStep\fR \fIt\fR \fItstep\fR \fIxvec\fR \fIfunc\fR Set a single step in the numerical integration of a system of differential equations\&. The method used is Euler's\&. .RS .TP \fIt\fR Value of the independent variable (typically time) at the beginning of the step\&. .TP \fItstep\fR Step size for the independent variable\&. .TP \fIxvec\fR List (vector) of dependent values .TP \fIfunc\fR 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)\&. .RE .sp .TP \fB::math::calculus::heunStep\fR \fIt\fR \fItstep\fR \fIxvec\fR \fIfunc\fR Set a single step in the numerical integration of a system of differential equations\&. The method used is Heun's\&. .RS .TP \fIt\fR Value of the independent variable (typically time) at the beginning of the step\&. .TP \fItstep\fR Step size for the independent variable\&. .TP \fIxvec\fR List (vector) of dependent values .TP \fIfunc\fR 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)\&. .RE .sp .TP \fB::math::calculus::rungeKuttaStep\fR \fIt\fR \fItstep\fR \fIxvec\fR \fIfunc\fR Set a single step in the numerical integration of a system of differential equations\&. The method used is Runge-Kutta 4th order\&. .RS .TP \fIt\fR Value of the independent variable (typically time) at the beginning of the step\&. .TP \fItstep\fR Step size for the independent variable\&. .TP \fIxvec\fR List (vector) of dependent values .TP \fIfunc\fR 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)\&. .RE .sp .TP \fB::math::calculus::boundaryValueSecondOrder\fR \fIcoeff_func\fR \fIforce_func\fR \fIleftbnd\fR \fIrightbnd\fR \fInostep\fR Solve a second order linear differential equation with boundary values at two sides\&. The equation has to be of the form (the "conservative" form): .CS d dy d -- A(x)-- + -- B(x)y + C(x)y = D(x) dx dx dx .CE .IP Ordinarily, such an equation would be written as: .CS d2y dy a(x)--- + b(x)-- + c(x) y = D(x) dx2 dx .CE .IP 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: .CS A(x) = a(x) B(x) = b(x) - a'(x) C(x) = c(x) - B'(x) = c(x) - b'(x) + a''(x) .CE .IP Because of the differentiation, however, it is much easier to ask the user to provide the functions A, B and C directly\&. .RS .TP \fIcoeff_func\fR Procedure returning the three coefficients (A, B, C) of the equation, taking as its one argument the x-coordinate\&. .TP \fIforce_func\fR Procedure returning the right-hand side (D) as a function of the x-coordinate\&. .TP \fIleftbnd\fR A list of two values: the x-coordinate of the left boundary and the value at that boundary\&. .TP \fIrightbnd\fR A list of two values: the x-coordinate of the right boundary and the value at that boundary\&. .TP \fInostep\fR Number of steps by which to discretise the interval\&. The procedure returns a list of x-coordinates and the approximated values of the solution\&. .RE .sp .TP \fB::math::calculus::solveTriDiagonal\fR \fIacoeff\fR \fIbcoeff\fR \fIccoeff\fR \fIdvalue\fR Solve a system of linear equations Ax = b with A a tridiagonal matrix\&. Returns the solution as a list\&. .RS .TP \fIacoeff\fR List of values on the lower diagonal .TP \fIbcoeff\fR List of values on the main diagonal .TP \fIccoeff\fR List of values on the upper diagonal .TP \fIdvalue\fR List of values on the righthand-side .RE .sp .TP \fB::math::calculus::newtonRaphson\fR \fIfunc\fR \fIderiv\fR \fIinitval\fR Determine the root of an equation given by .CS func(x) = 0 .CE .IP using the method of Newton-Raphson\&. The procedure takes the following arguments: .RS .TP \fIfunc\fR Procedure that returns the value the function at x .TP \fIderiv\fR Procedure that returns the derivative of the function at x .TP \fIinitval\fR Initial value for x .RE .sp .TP \fB::math::calculus::newtonRaphsonParameters\fR \fImaxiter\fR \fItolerance\fR Set the numerical parameters for the Newton-Raphson method: .RS .TP \fImaxiter\fR Maximum number of iteration steps (defaults to 20) .TP \fItolerance\fR Relative precision (defaults to 0\&.001) .RE .TP \fB::math::calculus::regula_falsi\fR \fIf\fR \fIxb\fR \fIxe\fR \fIeps\fR 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\&. .sp The method used is the so-called \fIregula falsi\fR or \fIfalse position\fR 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\&. .sp In contrast to Newton-Raphson there is no need for the computation of the function's derivative\&. .RS .TP command \fIf\fR Name of the command that evaluates the function for which the zero is to be returned .TP float \fIxb\fR Start of the interval in which the zero is supposed to lie .TP float \fIxe\fR End of the interval .TP float \fIeps\fR Relative allowed error (defaults to 1\&.0e-4) .RE .PP .PP \fINotes:\fR .PP Several of the above procedures take the \fInames\fR of procedures as arguments\&. To avoid problems with the \fIvisibility\fR 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: .CS 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] } } .CE .PP Enhancements for the second-order boundary value problem: .IP \(bu Other types of boundary conditions (zero gradient, zero flux) .IP \(bu Other schematisation of the first-order term (now central differences are used, but upstream differences might be useful too)\&. .PP .SH EXAMPLES Let us take a few simple examples: .PP Integrate x over the interval [0,100] (20 steps): .CS proc linear_func { x } { return $x } puts "Integral: [::math::calculus::integral 0 100 20 linear_func]" .CE For simple functions, the alternative could be: .CS puts "Integral: [::math::calculus::integralExpr 0 100 20 {$x}]" .CE Do not forget the braces! .PP The differential equation for a dampened oscillator: .PP .CS x'' + rx' + wx = 0 .CE .PP can be split into a system of first-order equations: .PP .CS x' = y y' = -ry - wx .CE .PP Then this system can be solved with code like this: .PP .CS 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 } .CE .PP Suppose we have the boundary value problem: .PP .CS Dy'' + ky = 0 x = 0: y = 1 x = L: y = 0 .CE .PP This boundary value problem could originate from the diffusion of a decaying substance\&. .PP It can be solved with the following fragment: .PP .CS 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] .CE .SH "BUGS, IDEAS, FEEDBACK" This document, and the package it describes, will undoubtedly contain bugs and other problems\&. Please report such in the category \fImath :: calculus\fR of the \fITcllib Trackers\fR [http://core\&.tcl\&.tk/tcllib/reportlist]\&. Please also report any ideas for enhancements you may have for either package and/or documentation\&. .SH "SEE ALSO" romberg .SH KEYWORDS calculus, differential equations, integration, math, roots .SH CATEGORY Mathematics .SH COPYRIGHT .nf Copyright (c) 2002,2003,2004 Arjen Markus .fi