/usr/share/genius/gel/calculus/integration.gel is in genius-common 1.0.23-3.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 | # Numerical integration
#
# The algorithms are described in:
# Numerical Analysis, 5th edition
# by Richard L. Burden and J. Douglas Faires
# PWS Publishing Company, Boston, 1993.
# Library of congress: QA 297 B84 1993
# In the below, f indicates the function whose integral we wish to determine,
# a,b indicate the left and right endpoints of the interval over which
# we wish to integrate, and n is the number of intervals into which we
# divide [a,b]
# These methods all return one value, the value of the integral
# Currently only works for real functions of a real variable
# Composite Simpson's Rule, Section 4.4, Algorithm 4.1, p. 186
# Note that this has error term = max(f'''')*h^4*(b-a)/180,
# where h=(b-a)/n
# If we can get maximums and derivatives, this would allow us to determine
# automatically what n should be.
# Composite simpson's rule is implemented as a built in function
SetHelp ("CompositeSimpsonsRuleTolerance", "calculus", "Integration of f by Composite Simpson's Rule on the interval [a,b] with the number of steps calculated by the fourth derivative bound and the desired tolerance")
function CompositeSimpsonsRuleTolerance(f,a,b,FourthDerivativeBound,Tolerance) =
(
local *;
# Error term = max(f'''')*h^4*(b-a)/180,
# where h=(b-a)/n
n = ceil(|FourthDerivativeBound*(b-a)^5 / (180*Tolerance)|^(1/4));
# Note that this is done automatically by CompositeSimpsonsRule
# function:
#if IsOdd(n) then n = n+1;
CompositeSimpsonsRule (f, a, b, n)
)
# FIXME: reference, though this is really stupid
SetHelp ("MidpointRule", "calculus", "Integration by midpoint rule")
function MidpointRule(f,a,b,n) =
(
local *;
if(not IsFunction(f)) then
(error("MidpointRule: argument 1 must be a function");bailout)
else if(not IsReal(a) or not IsReal(b)) then
(error("MidpointRule: arguments 2, 3 must be real values");bailout)
else if(not IsInteger(n)) then
(error("MidpointRule: argument 4 must be an integer");bailout);
## check bounds
if(a>b) then (error("MidpointRule: argument 2 must be less than or equal to argument 3");bailout)
else if(n<= 0) then (error("MidpointRule: argument 4 must be positive");bailout);
len = float(b-a);
s = sum i=1 to n do float(f(a+(len*(i-0.5))/n));
(s*len)/n
)
SetHelp ("NumericalIntegralSteps", "parameters", "Steps to perform in NumericalIntegral")
parameter NumericalIntegralSteps = 1000
SetHelp ("NumericalIntegralFunction", "parameters", "The function used for numerical integration in NumericalIntegral")
parameter NumericalIntegralFunction = `CompositeSimpsonsRule
SetHelp ("NumericalIntegral", "calculus", "Integration by rule set in NumericalIntegralFunction of f from a to b using NumericalIntegralSteps steps")
function NumericalIntegral(f,a,b) = (
local *;
NumericalIntegralFunction call (f,a,b,NumericalIntegralSteps)
)
|