Numerical integration: Difference between revisions

Content added Content deleted
m (syntax highlighting fixup automation)
Line 48: Line 48:
{{trans|Nim}}
{{trans|Nim}}


<lang 11l>F left_rect((Float -> Float) f, Float x, Float h) -> Float
<syntaxhighlight lang="11l">F left_rect((Float -> Float) f, Float x, Float h) -> Float
R f(x)
R f(x)


Line 87: Line 87:
(simpson, ‘simpson’)]
(simpson, ‘simpson’)]
print("#. integrated using #.\n from #. to #. (#. steps) = #.".format(
print("#. integrated using #.\n from #. to #. (#. steps) = #.".format(
func_name, rule_name, a, b, steps, integrate(func, a, b, steps, rule)))</lang>
func_name, rule_name, a, b, steps, integrate(func, a, b, steps, rule)))</syntaxhighlight>


{{out}}
{{out}}
Line 135: Line 135:
=={{header|ActionScript}}==
=={{header|ActionScript}}==
Integration functions:
Integration functions:
<lang ActionScript>function leftRect(f:Function, a:Number, b:Number, n:uint):Number
<syntaxhighlight lang="actionscript">function leftRect(f:Function, a:Number, b:Number, n:uint):Number
{
{
var sum:Number = 0;
var sum:Number = 0;
Line 185: Line 185:
}
}
return (dx/6) * (f(a) + f(b) + 4*sum1 + 2*sum2);
return (dx/6) * (f(a) + f(b) + 4*sum1 + 2*sum2);
}</lang>
}</syntaxhighlight>
Usage:
Usage:
<lang ActionScript>function f1(n:Number):Number {
<syntaxhighlight lang="actionscript">function f1(n:Number):Number {
return (2/(1+ 4*(n*n)));
return (2/(1+ 4*(n*n)));
}
}
Line 195: Line 195:
trace(trapezium(f1, -1, 2 ,4 ));
trace(trapezium(f1, -1, 2 ,4 ));
trace(simpson(f1, -1, 2 ,4 ));
trace(simpson(f1, -1, 2 ,4 ));
</syntaxhighlight>
</lang>


=={{header|Ada}}==
=={{header|Ada}}==
Specification of a generic package implementing the five specified kinds of numerical integration:
Specification of a generic package implementing the five specified kinds of numerical integration:
<lang ada>generic
<syntaxhighlight lang="ada">generic
type Scalar is digits <>;
type Scalar is digits <>;
with function F (X : Scalar) return Scalar;
with function F (X : Scalar) return Scalar;
Line 208: Line 208:
function Trapezium (A, B : Scalar; N : Positive) return Scalar;
function Trapezium (A, B : Scalar; N : Positive) return Scalar;
function Simpsons (A, B : Scalar; N : Positive) return Scalar;
function Simpsons (A, B : Scalar; N : Positive) return Scalar;
end Integrate;</lang>
end Integrate;</syntaxhighlight>
An alternative solution is to pass a function reference to the integration function. This solution is probably slightly faster, and works even with Ada83. One could also make each integration function generic, instead of making the whole package generic.
An alternative solution is to pass a function reference to the integration function. This solution is probably slightly faster, and works even with Ada83. One could also make each integration function generic, instead of making the whole package generic.


Body of the package implementing numerical integration:
Body of the package implementing numerical integration:
<lang ada>package body Integrate is
<syntaxhighlight lang="ada">package body Integrate is
function Left_Rectangular (A, B : Scalar; N : Positive) return Scalar is
function Left_Rectangular (A, B : Scalar; N : Positive) return Scalar is
H : constant Scalar := (B - A) / Scalar (N);
H : constant Scalar := (B - A) / Scalar (N);
Line 275: Line 275:
return (H / 3.0) * (F (A) + F (B) + 4.0 * Sum_U + 2.0 * Sum_E);
return (H / 3.0) * (F (A) + F (B) + 4.0 * Sum_U + 2.0 * Sum_E);
end Simpsons;
end Simpsons;
end Integrate;</lang>
end Integrate;</syntaxhighlight>


Test driver:
Test driver:
<lang ada>with Ada.Text_IO, Ada.Integer_Text_IO;
<syntaxhighlight lang="ada">with Ada.Text_IO, Ada.Integer_Text_IO;
with Integrate;
with Integrate;


Line 382: Line 382:
end X;
end X;
end Numerical_Integration;
end Numerical_Integration;
</syntaxhighlight>
</lang>


=={{header|ALGOL 68}}==
=={{header|ALGOL 68}}==
<lang algol68>MODE F = PROC(LONG REAL)LONG REAL;
<syntaxhighlight lang="algol68">MODE F = PROC(LONG REAL)LONG REAL;


###############
###############
Line 502: Line 502:
test integrators( "x ", ( LONG REAL x )LONG REAL: x, 0, 6 000, 6 000 000 );
test integrators( "x ", ( LONG REAL x )LONG REAL: x, 0, 6 000, 6 000 000 );


SKIP</lang>
SKIP</syntaxhighlight>
{{out}}
{{out}}
<pre>
<pre>
Line 513: Line 513:
=={{header|ALGOL W}}==
=={{header|ALGOL W}}==
{{Trans|ALGOL 68}}
{{Trans|ALGOL 68}}
<lang algolw>begin % compare some numeric integration methods %
<syntaxhighlight lang="algolw">begin % compare some numeric integration methods %


long real procedure leftRect ( long real procedure f
long real procedure leftRect ( long real procedure f
Line 641: Line 641:
testIntegrators2( "x ", xValue, 0, 6000, 6000000 )
testIntegrators2( "x ", xValue, 0, 6000, 6000000 )
end
end
end.</lang>
end.</syntaxhighlight>


=={{header|AutoHotkey}}==
=={{header|AutoHotkey}}==
ahk [http://www.autohotkey.com/forum/viewtopic.php?t=44657&postdays=0&postorder=asc&start=139 discussion]
ahk [http://www.autohotkey.com/forum/viewtopic.php?t=44657&postdays=0&postorder=asc&start=139 discussion]
<lang autohotkey>MsgBox % Rect("fun", 0, 1, 10,-1) ; 0.45 left
<syntaxhighlight lang="autohotkey">MsgBox % Rect("fun", 0, 1, 10,-1) ; 0.45 left
MsgBox % Rect("fun", 0, 1, 10) ; 0.50 mid
MsgBox % Rect("fun", 0, 1, 10) ; 0.50 mid
MsgBox % Rect("fun", 0, 1, 10, 1) ; 0.55 right
MsgBox % Rect("fun", 0, 1, 10, 1) ; 0.55 right
Line 679: Line 679:
fun(x) { ; linear test function
fun(x) { ; linear test function
Return x
Return x
}</lang>
}</syntaxhighlight>


=={{header|BASIC}}==
=={{header|BASIC}}==
{{works with|QuickBasic|4.5}}
{{works with|QuickBasic|4.5}}
{{trans|Java}}
{{trans|Java}}
<lang qbasic>FUNCTION leftRect(a, b, n)
<syntaxhighlight lang="qbasic">FUNCTION leftRect(a, b, n)
h = (b - a) / n
h = (b - a) / n
sum = 0
sum = 0
Line 734: Line 734:


simpson = h / 6 * (f(a) + f(b) + 4 * sum1 + 2 * sum2)
simpson = h / 6 * (f(a) + f(b) + 4 * sum1 + 2 * sum2)
END FUNCTION</lang>
END FUNCTION</syntaxhighlight>


=={{header|BBC BASIC}}==
=={{header|BBC BASIC}}==
<lang bbcbasic> *FLOAT64
<syntaxhighlight lang="bbcbasic"> *FLOAT64
@% = 12 : REM Column width
@% = 12 : REM Column width
Line 807: Line 807:
NEXT
NEXT
x = a
x = a
= (d / 6) * (f + EVAL(x$) + 4 * s1 + 2 * s2)</lang>
= (d / 6) * (f + EVAL(x$) + 4 * s1 + 2 * s2)</syntaxhighlight>
'''Output:'''
'''Output:'''
<pre>
<pre>
Line 818: Line 818:


=={{header|C}}==
=={{header|C}}==
<lang c>#include <stdio.h>
<syntaxhighlight lang="c">#include <stdio.h>
#include <stdlib.h>
#include <stdlib.h>
#include <math.h>
#include <math.h>
Line 875: Line 875:


return h / 6.0 * (func(from) + func(to) + 4.0 * sum1 + 2.0 * sum2);
return h / 6.0 * (func(from) + func(to) + 4.0 * sum1 + 2.0 * sum2);
}</lang>
}</syntaxhighlight>


<lang c>/* test */
<syntaxhighlight lang="c">/* test */
double f3(double x)
double f3(double x)
{
{
Line 947: Line 947:
printf("\n");
printf("\n");
}
}
}</lang>
}</syntaxhighlight>


=={{header|C sharp|C#}}==
=={{header|C sharp|C#}}==
<lang csharp>using System;
<syntaxhighlight lang="csharp">using System;
using System.Collections.Generic;
using System.Collections.Generic;
using System.Linq;
using System.Linq;
Line 1,068: Line 1,068:
TestApproximationMethods(new DefiniteIntegral(x => x, new Interval(0, 6000)), 6000000);
TestApproximationMethods(new DefiniteIntegral(x => x, new Interval(0, 6000)), 6000000);
}
}
}</lang>
}</syntaxhighlight>
Output:
Output:
<lang>0.2499500025
<syntaxhighlight lang="text">0.2499500025
0.24999999875
0.24999999875
0.2500500025
0.2500500025
Line 1,089: Line 1,089:
18000003
18000003
18000000
18000000
18000000</lang>
18000000</syntaxhighlight>


=={{header|C++}}==
=={{header|C++}}==


Due to their similarity, it makes sense to make the integration method a policy.
Due to their similarity, it makes sense to make the integration method a policy.
<lang cpp>// the integration routine
<syntaxhighlight lang="cpp">// the integration routine
template<typename Method, typename F, typename Float>
template<typename Method, typename F, typename Float>
double integrate(F f, Float a, Float b, int steps, Method m)
double integrate(F f, Float a, Float b, int steps, Method m)
Line 1,156: Line 1,156:
double rr = integrate(f, 0.0, 1.0, 10, rectangular(rectangular::right));
double rr = integrate(f, 0.0, 1.0, 10, rectangular(rectangular::right));
double t = integrate(f, 0.0, 1.0, 10, trapezium());
double t = integrate(f, 0.0, 1.0, 10, trapezium());
double s = integrate(f, 0.0, 1.0, 10, simpson());</lang>
double s = integrate(f, 0.0, 1.0, 10, simpson());</syntaxhighlight>


=={{header|Chapel}}==
=={{header|Chapel}}==
<lang chapel>
<syntaxhighlight lang="chapel">
proc f1(x:real):real {
proc f1(x:real):real {
return x**3;
return x**3;
Line 1,289: Line 1,289:
writeln("simpsonsIntegration: calculated = ", calculated, "; exact = ", exact, "; difference = ", abs(calculated - exact));
writeln("simpsonsIntegration: calculated = ", calculated, "; exact = ", exact, "; difference = ", abs(calculated - exact));
writeln();
writeln();
</syntaxhighlight>
</lang>
output
output
<lang>
<syntaxhighlight lang="text">
f(x) = x**3 with 100 steps from 0 to 1
f(x) = x**3 with 100 steps from 0 to 1
leftRectangleIntegration: calculated = 0.245025; exact = 0.25; difference = 0.004975
leftRectangleIntegration: calculated = 0.245025; exact = 0.25; difference = 0.004975
Line 1,319: Line 1,319:
trapezoidIntegration: calculated = 1.8e+07; exact = 1.8e+07; difference = 3.72529e-09
trapezoidIntegration: calculated = 1.8e+07; exact = 1.8e+07; difference = 3.72529e-09
simpsonsIntegration: calculated = 1.8e+07; exact = 1.8e+07; difference = 0.0
simpsonsIntegration: calculated = 1.8e+07; exact = 1.8e+07; difference = 0.0
</syntaxhighlight>
</lang>


=={{header|CoffeeScript}}==
=={{header|CoffeeScript}}==
{{trans|python}}
{{trans|python}}
<lang coffeescript>
<syntaxhighlight lang="coffeescript">
rules =
rules =
left_rect: (f, x, h) -> f(x)
left_rect: (f, x, h) -> f(x)
Line 1,357: Line 1,357:
result = integrate func, a, b, steps, rule
result = integrate func, a, b, steps, rule
console.log rule_name, result
console.log rule_name, result
</syntaxhighlight>
</lang>
output
output
<lang>
<syntaxhighlight lang="text">
> coffee numerical_integration.coffee
> coffee numerical_integration.coffee
-- tests for cube with 100 steps from 0 to 1
-- tests for cube with 100 steps from 0 to 1
Line 1,385: Line 1,385:
trapezium 17999999.999999993
trapezium 17999999.999999993
simpson 17999999.999999993
simpson 17999999.999999993
</syntaxhighlight>
</lang>


=={{header|Comal}}==
=={{header|Comal}}==
{{works with|OpenComal on Linux}}
{{works with|OpenComal on Linux}}
<syntaxhighlight lang="comal">
<lang Comal>
1000 PRINT "F(X)";" FROM";" TO";" L-Rect";" M-Rect";" R-Rect ";" Trapez";" Simpson"
1000 PRINT "F(X)";" FROM";" TO";" L-Rect";" M-Rect";" R-Rect ";" Trapez";" Simpson"
1010 fromval:=0
1010 fromval:=0
Line 1,484: Line 1,484:
1920 RETURN x
1920 RETURN x
1930 ENDFUNC
1930 ENDFUNC
</syntaxhighlight>
</lang>
{{out}}
{{out}}
<pre>
<pre>
Line 1,496: Line 1,496:
=={{header|Common Lisp}}==
=={{header|Common Lisp}}==


<lang lisp>(defun left-rectangle (f a b n &aux (d (/ (- b a) n)))
<syntaxhighlight lang="lisp">(defun left-rectangle (f a b n &aux (d (/ (- b a) n)))
(* d (loop for x from a below b by d summing (funcall f x))))
(* d (loop for x from a below b by d summing (funcall f x))))


Line 1,522: Line 1,522:
(funcall f b)
(funcall f b)
(* 4 sum1)
(* 4 sum1)
(* 2 sum2))))))</lang>
(* 2 sum2))))))</syntaxhighlight>


=={{header|D}}==
=={{header|D}}==
<lang d>import std.stdio, std.typecons, std.typetuple;
<syntaxhighlight lang="d">import std.stdio, std.typecons, std.typetuple;


template integrate(alias method) {
template integrate(alias method) {
Line 1,587: Line 1,587:
writeln();
writeln();
}
}
}</lang>
}</syntaxhighlight>
Output:
Output:
<pre>rectangular left: 0.202500
<pre>rectangular left: 0.202500
Line 1,614: Line 1,614:
===A faster version===
===A faster version===
This version avoids function pointers and delegates, same output:
This version avoids function pointers and delegates, same output:
<lang d>import std.stdio, std.typecons, std.typetuple;
<syntaxhighlight lang="d">import std.stdio, std.typecons, std.typetuple;


template integrate(alias method) {
template integrate(alias method) {
Line 1,684: Line 1,684:
writeln();
writeln();
}
}
}</lang>
}</syntaxhighlight>
=={{header|Delphi}}==
=={{header|Delphi}}==
{{libheader| System.SysUtils}}
{{libheader| System.SysUtils}}
{{Trans|Python}}
{{Trans|Python}}
<lang Delphi>program Numerical_integration;
<syntaxhighlight lang="delphi">program Numerical_integration;


{$APPTYPE CONSOLE}
{$APPTYPE CONSOLE}
Line 1,785: Line 1,785:
end;
end;
readln;
readln;
end.</lang>
end.</syntaxhighlight>
{{out}}
{{out}}
<pre>Integrate x^3
<pre>Integrate x^3
Line 1,815: Line 1,815:
{{trans|Python}}
{{trans|Python}}


<lang e>pragma.enable("accumulator")
<syntaxhighlight lang="e">pragma.enable("accumulator")


def leftRect(f, x, h) {
def leftRect(f, x, h) {
Line 1,840: Line 1,840:
def h := (b-a) / steps
def h := (b-a) / steps
return h * accum 0 for i in 0..!steps { _ + meth(f, a+i*h, h) }
return h * accum 0 for i in 0..!steps { _ + meth(f, a+i*h, h) }
}</lang>
}</syntaxhighlight>
<lang e>? integrate(fn x { x ** 2 }, 3.0, 7.0, 30, simpson)
<syntaxhighlight lang="e">? integrate(fn x { x ** 2 }, 3.0, 7.0, 30, simpson)
# value: 105.33333333333334
# value: 105.33333333333334


? integrate(fn x { x ** 9 }, 0, 1, 300, simpson)
? integrate(fn x { x ** 9 }, 0, 1, 300, simpson)
# value: 0.10000000002160479</lang>
# value: 0.10000000002160479</syntaxhighlight>


=={{header|Elixir}}==
=={{header|Elixir}}==
<lang elixir>defmodule Numerical do
<syntaxhighlight lang="elixir">defmodule Numerical do
@funs ~w(leftrect midrect rightrect trapezium simpson)a
@funs ~w(leftrect midrect rightrect trapezium simpson)a
Line 1,884: Line 1,884:
f4 = fn x -> x end
f4 = fn x -> x end
IO.puts "\nf(x) = x, where x is [0,6000], with 6,000,000 approximations."
IO.puts "\nf(x) = x, where x is [0,6000], with 6,000,000 approximations."
Numerical.integrate(f4, 0, 6000, 6_000_000)</lang>
Numerical.integrate(f4, 0, 6000, 6_000_000)</syntaxhighlight>


{{out}}
{{out}}
Line 1,918: Line 1,918:


=={{header|Euphoria}}==
=={{header|Euphoria}}==
<lang euphoria>function int_leftrect(sequence bounds, integer n, integer func_id)
<syntaxhighlight lang="euphoria">function int_leftrect(sequence bounds, integer n, integer func_id)
atom h, sum
atom h, sum
h = (bounds[2]-bounds[1])/n
h = (bounds[2]-bounds[1])/n
Line 1,996: Line 1,996:
? int_rightrect({0,10},1000,routine_id("x"))
? int_rightrect({0,10},1000,routine_id("x"))
? int_midrect({0,10},1000,routine_id("x"))
? int_midrect({0,10},1000,routine_id("x"))
? int_simpson({0,10},1000,routine_id("x"))</lang>
? int_simpson({0,10},1000,routine_id("x"))</syntaxhighlight>


Output:
Output:
Line 2,016: Line 2,016:


=={{header|F Sharp}}==
=={{header|F Sharp}}==
<lang fsharp>
<syntaxhighlight lang="fsharp">
// integration methods
// integration methods
let left f dx x = f x * dx
let left f dx x = f x * dx
Line 2,042: Line 2,042:
|> Seq.map (fun ((f, a, b, n), method) -> integrate a b f n method)
|> Seq.map (fun ((f, a, b, n), method) -> integrate a b f n method)
|> Seq.iter (printfn "%f")
|> Seq.iter (printfn "%f")
</syntaxhighlight>
</lang>


=={{header|Factor}}==
=={{header|Factor}}==
<lang factor>
<syntaxhighlight lang="factor">
USE: math.functions
USE: math.functions
IN: scratchpad 0 1 [ 3 ^ ] integrate-simpson .
IN: scratchpad 0 1 [ 3 ^ ] integrate-simpson .
Line 2,057: Line 2,057:
IN: scratchpad 6000000 num-steps set-global
IN: scratchpad 6000000 num-steps set-global
IN: scratchpad 0 6000 [ ] integrate-simpson .
IN: scratchpad 0 6000 [ ] integrate-simpson .
18000000</lang>
18000000</syntaxhighlight>


=={{header|Forth}}==
=={{header|Forth}}==
<lang forth>fvariable step
<syntaxhighlight lang="forth">fvariable step


defer method ( fn F: x -- fn[x] )
defer method ( fn F: x -- fn[x] )
Line 2,099: Line 2,099:
test mid fn2 \ 2.496091
test mid fn2 \ 2.496091
test trap fn2 \ 2.351014
test trap fn2 \ 2.351014
test simpson fn2 \ 2.447732</lang>
test simpson fn2 \ 2.447732</syntaxhighlight>


=={{header|Fortran}}==
=={{header|Fortran}}==
In ISO Fortran 95 and later if function f() is not already defined to be "elemental", define an elemental wrapper function around it to allow for array-based initialization:
In ISO Fortran 95 and later if function f() is not already defined to be "elemental", define an elemental wrapper function around it to allow for array-based initialization:
<lang fortran>elemental function elemf(x)
<syntaxhighlight lang="fortran">elemental function elemf(x)
real :: elemf, x
real :: elemf, x
elemf = f(x)
elemf = f(x)
end function elemf</lang>
end function elemf</syntaxhighlight>


Use Array Initializers, Pointers, Array invocation of Elemental functions, Elemental array-array and array-scalar arithmetic, and the SUM intrinsic function. Methods are collected into a single function in a module.
Use Array Initializers, Pointers, Array invocation of Elemental functions, Elemental array-array and array-scalar arithmetic, and the SUM intrinsic function. Methods are collected into a single function in a module.
<lang fortran>module Integration
<syntaxhighlight lang="fortran">module Integration
implicit none
implicit none


Line 2,183: Line 2,183:
end function integrate
end function integrate


end module Integration</lang>
end module Integration</syntaxhighlight>


Usage example:
Usage example:
<lang fortran>program IntegrationTest
<syntaxhighlight lang="fortran">program IntegrationTest
use Integration
use Integration
use FunctionHolder
use FunctionHolder
Line 2,197: Line 2,197:
print *, integrate(afun, 0., 3**(1/3.), method='trapezoid')
print *, integrate(afun, 0., 3**(1/3.), method='trapezoid')


end program IntegrationTest</lang>
end program IntegrationTest</syntaxhighlight>


The FunctionHolder module:
The FunctionHolder module:


<lang fortran>module FunctionHolder
<syntaxhighlight lang="fortran">module FunctionHolder
implicit none
implicit none
Line 2,213: Line 2,213:
end function afun
end function afun
end module FunctionHolder</lang>
end module FunctionHolder</syntaxhighlight>


=={{header|FreeBASIC}}==
=={{header|FreeBASIC}}==
Based on the BASIC entry and the BBC BASIC entry
Based on the BASIC entry and the BBC BASIC entry
<lang freebasic>' version 17-09-2015
<syntaxhighlight lang="freebasic">' version 17-09-2015
' compile with: fbc -s console
' compile with: fbc -s console


Line 2,341: Line 2,341:
Print : Print "hit any key to end program"
Print : Print "hit any key to end program"
Sleep
Sleep
End</lang>
End</syntaxhighlight>
{{out}}
{{out}}
<pre>function range steps leftrect midrect rightrect trap simpson
<pre>function range steps leftrect midrect rightrect trap simpson
Line 2,350: Line 2,350:


=={{header|Go}}==
=={{header|Go}}==
<lang go>package main
<syntaxhighlight lang="go">package main


import (
import (
Line 2,513: Line 2,513:
fmt.Println("")
fmt.Println("")
}
}
}</lang>
}</syntaxhighlight>
{{out}}
{{out}}
<pre>
<pre>
Line 2,554: Line 2,554:
=={{header|Groovy}}==
=={{header|Groovy}}==
Solution:
Solution:
<lang groovy>def assertBounds = { List bounds, int nRect ->
<syntaxhighlight lang="groovy">def assertBounds = { List bounds, int nRect ->
assert (bounds.size() == 2) && (bounds[0] instanceof Double) && (bounds[1] instanceof Double) && (nRect > 0)
assert (bounds.size() == 2) && (bounds[0] instanceof Double) && (bounds[1] instanceof Double) && (nRect > 0)
}
}
Line 2,597: Line 2,597:
h/3*((fLeft + fRight).sum() + 4*(fMid.sum()))
h/3*((fLeft + fRight).sum() + 4*(fMid.sum()))
}
}
}</lang>
}</syntaxhighlight>


Test:
Test:


Each "nRect" (number of rectangles) value given below is the minimum value that meets the tolerance condition for the given circumstances (function-to-integrate, integral-type and integral-bounds).
Each "nRect" (number of rectangles) value given below is the minimum value that meets the tolerance condition for the given circumstances (function-to-integrate, integral-type and integral-bounds).
<lang groovy>double tolerance = 0.0001 // allowable "wrongness", ensures accuracy to 1 in 10,000
<syntaxhighlight lang="groovy">double tolerance = 0.0001 // allowable "wrongness", ensures accuracy to 1 in 10,000


double sinIntegralCalculated = -(Math.cos(Math.PI) - Math.cos(0d))
double sinIntegralCalculated = -(Math.cos(Math.PI) - Math.cos(0d))
Line 2,641: Line 2,641:
assert ((simpsonsIntegral([0d, Math.PI], 1, cubicPoly) - cpIntegralCalc0ToPI)/ cpIntegralCalc0ToPI).abs() < tolerance**2.75 // 1 in 100 billion
assert ((simpsonsIntegral([0d, Math.PI], 1, cubicPoly) - cpIntegralCalc0ToPI)/ cpIntegralCalc0ToPI).abs() < tolerance**2.75 // 1 in 100 billion
double cpIntegralCalcMinusEToPI = (cubicPolyAntiDeriv(Math.PI) - cubicPolyAntiDeriv(-Math.E))
double cpIntegralCalcMinusEToPI = (cubicPolyAntiDeriv(Math.PI) - cubicPolyAntiDeriv(-Math.E))
assert ((simpsonsIntegral([-Math.E, Math.PI], 1, cubicPoly) - cpIntegralCalcMinusEToPI)/ cpIntegralCalcMinusEToPI).abs() < tolerance**2.5 // 1 in 10 billion</lang>
assert ((simpsonsIntegral([-Math.E, Math.PI], 1, cubicPoly) - cpIntegralCalcMinusEToPI)/ cpIntegralCalcMinusEToPI).abs() < tolerance**2.5 // 1 in 10 billion</syntaxhighlight>


Requested Demonstrations:
Requested Demonstrations:
<lang groovy>println "f(x) = x**3, where x is [0,1], with 100 approximations. The exact result is 1/4, or 0.25."
<syntaxhighlight lang="groovy">println "f(x) = x**3, where x is [0,1], with 100 approximations. The exact result is 1/4, or 0.25."
println ([" LeftRect": leftRectIntegral([0d, 1d], 100) { it**3 }])
println ([" LeftRect": leftRectIntegral([0d, 1d], 100) { it**3 }])
println (["RightRect": rightRectIntegral([0d, 1d], 100) { it**3 }])
println (["RightRect": rightRectIntegral([0d, 1d], 100) { it**3 }])
Line 2,674: Line 2,674:
println (["Trapezoid": trapezoidIntegral([0d, 6000d], 6000000) { it }])
println (["Trapezoid": trapezoidIntegral([0d, 6000d], 6000000) { it }])
println ([" Simpsons": simpsonsIntegral([0d, 6000d], 6000000) { it }])
println ([" Simpsons": simpsonsIntegral([0d, 6000d], 6000000) { it }])
println ()</lang>
println ()</syntaxhighlight>


Output:
Output:
Line 2,709: Line 2,709:
Different approach from most of the other examples: First, the function ''f'' might be expensive to calculate, and so it should not be evaluated several times. So, ideally, we want to have positions ''x'' and weights ''w'' for each method and then just calculate the approximation of the integral by
Different approach from most of the other examples: First, the function ''f'' might be expensive to calculate, and so it should not be evaluated several times. So, ideally, we want to have positions ''x'' and weights ''w'' for each method and then just calculate the approximation of the integral by


<lang haskell>approx f xs ws = sum [w * f x | (x,w) <- zip xs ws]</lang>
<syntaxhighlight lang="haskell">approx f xs ws = sum [w * f x | (x,w) <- zip xs ws]</syntaxhighlight>


Second, let's to generalize all integration methods into one scheme. The methods can all be characterized by the coefficients ''vs'' they use in a particular interval. These will be fractions, and for terseness, we extract the denominator as an extra argument ''v''.
Second, let's to generalize all integration methods into one scheme. The methods can all be characterized by the coefficients ''vs'' they use in a particular interval. These will be fractions, and for terseness, we extract the denominator as an extra argument ''v''.
Line 2,715: Line 2,715:
Now there are the closed formulas (which include the endpoints) and the open formulas (which exclude them). Let's do the open formulas first, because then the coefficients don't overlap:
Now there are the closed formulas (which include the endpoints) and the open formulas (which exclude them). Let's do the open formulas first, because then the coefficients don't overlap:
<lang haskell>integrateOpen :: Fractional a => a -> [a] -> (a -> a) -> a -> a -> Int -> a
<syntaxhighlight lang="haskell">integrateOpen :: Fractional a => a -> [a] -> (a -> a) -> a -> a -> Int -> a
integrateOpen v vs f a b n = approx f xs ws * h / v where
integrateOpen v vs f a b n = approx f xs ws * h / v where
m = fromIntegral (length vs) * n
m = fromIntegral (length vs) * n
Line 2,721: Line 2,721:
ws = concat $ replicate n vs
ws = concat $ replicate n vs
c = a + h/2
c = a + h/2
xs = [c + h * fromIntegral i | i <- [0..m-1]]</lang>
xs = [c + h * fromIntegral i | i <- [0..m-1]]</syntaxhighlight>


Similarly for the closed formulas, but we need an additional function ''overlap'' which sums the coefficients overlapping at the interior interval boundaries:
Similarly for the closed formulas, but we need an additional function ''overlap'' which sums the coefficients overlapping at the interior interval boundaries:
<lang haskell>integrateClosed :: Fractional a => a -> [a] -> (a -> a) -> a -> a -> Int -> a
<syntaxhighlight lang="haskell">integrateClosed :: Fractional a => a -> [a] -> (a -> a) -> a -> a -> Int -> a
integrateClosed v vs f a b n = approx f xs ws * h / v where
integrateClosed v vs f a b n = approx f xs ws * h / v where
m = fromIntegral (length vs - 1) * n
m = fromIntegral (length vs - 1) * n
Line 2,738: Line 2,738:
inter n [] = x : inter (n-1) xs
inter n [] = x : inter (n-1) xs
inter n [y] = (x+y) : inter (n-1) xs
inter n [y] = (x+y) : inter (n-1) xs
inter n (y:ys) = y : inter n ys</lang>
inter n (y:ys) = y : inter n ys</syntaxhighlight>


And now we can just define
And now we can just define


<lang haskell>intLeftRect = integrateClosed 1 [1,0]
<syntaxhighlight lang="haskell">intLeftRect = integrateClosed 1 [1,0]
intRightRect = integrateClosed 1 [0,1]
intRightRect = integrateClosed 1 [0,1]
intMidRect = integrateOpen 1 [1]
intMidRect = integrateOpen 1 [1]
intTrapezium = integrateClosed 2 [1,1]
intTrapezium = integrateClosed 2 [1,1]
intSimpson = integrateClosed 3 [1,4,1]</lang>
intSimpson = integrateClosed 3 [1,4,1]</syntaxhighlight>


or, as easily, some additional schemes:
or, as easily, some additional schemes:


<lang haskell>intMilne = integrateClosed 45 [14,64,24,64,14]
<syntaxhighlight lang="haskell">intMilne = integrateClosed 45 [14,64,24,64,14]
intOpen1 = integrateOpen 2 [3,3]
intOpen1 = integrateOpen 2 [3,3]
intOpen2 = integrateOpen 3 [8,-4,8]</lang>
intOpen2 = integrateOpen 3 [8,-4,8]</syntaxhighlight>


Some examples:
Some examples:
Line 2,769: Line 2,769:
The whole program:
The whole program:


<lang haskell>approx
<syntaxhighlight lang="haskell">approx
:: Fractional a
:: Fractional a
=> (a1 -> a) -> [a1] -> [a] -> a
=> (a1 -> a) -> [a1] -> [a] -> a
Line 2,857: Line 2,857:
integrations
integrations
where
where
indent n = take n . (++ replicate n ' ')</lang>
indent n = take n . (++ replicate n ' ')</syntaxhighlight>
{{Out}}
{{Out}}
<pre>f(x) = x^3 [0.0,1.0] (100 approximations)
<pre>f(x) = x^3 [0.0,1.0] (100 approximations)
Line 2,890: Line 2,890:
=={{header|J}}==
=={{header|J}}==
===Solution:===
===Solution:===
<lang j>integrate=: adverb define
<syntaxhighlight lang="j">integrate=: adverb define
'a b steps'=. 3{.y,128
'a b steps'=. 3{.y,128
size=. (b - a)%steps
size=. (b - a)%steps
Line 2,900: Line 2,900:
trapezium=: adverb def '-: +/ u y'
trapezium=: adverb def '-: +/ u y'


simpson =: adverb def '6 %~ +/ 1 1 4 * u y, -:+/y'</lang>
simpson =: adverb def '6 %~ +/ 1 1 4 * u y, -:+/y'</syntaxhighlight>
===Example usage===
===Example usage===
====Required Examples====
====Required Examples====
<lang j> Ir=: rectangle integrate
<syntaxhighlight lang="j"> Ir=: rectangle integrate
It=: trapezium integrate
It=: trapezium integrate
Is=: simpson integrate
Is=: simpson integrate
Line 2,930: Line 2,930:
1.8e7
1.8e7
] Is 0 6000 6e6
] Is 0 6000 6e6
1.8e7</lang>
1.8e7</syntaxhighlight>
====Older Examples====
====Older Examples====
Integrate <code>square</code> (<code>*:</code>) from 0 to &pi; in 10 steps using various methods.
Integrate <code>square</code> (<code>*:</code>) from 0 to &pi; in 10 steps using various methods.
<lang j> *: rectangle integrate 0 1p1 10
<syntaxhighlight lang="j"> *: rectangle integrate 0 1p1 10
10.3095869962
10.3095869962
*: trapezium integrate 0 1p1 10
*: trapezium integrate 0 1p1 10
10.3871026879
10.3871026879
*: simpson integrate 0 1p1 10
*: simpson integrate 0 1p1 10
10.3354255601</lang>
10.3354255601</syntaxhighlight>
Integrate <code>sin</code> from 0 to &pi; in 10 steps using various methods.
Integrate <code>sin</code> from 0 to &pi; in 10 steps using various methods.
<lang j> sin=: 1&o.
<syntaxhighlight lang="j"> sin=: 1&o.
sin rectangle integrate 0 1p1 10
sin rectangle integrate 0 1p1 10
2.00824840791
2.00824840791
Line 2,946: Line 2,946:
1.98352353751
1.98352353751
sin simpson integrate 0 1p1 10
sin simpson integrate 0 1p1 10
2.00000678444</lang>
2.00000678444</syntaxhighlight>
===Aside===
===Aside===
Note that J has a primitive verb <code>p..</code> for integrating polynomials. For example the integral of <math>x^2</math> (which can be described in terms of its coefficients as <code>0 0 1</code>) is:
Note that J has a primitive verb <code>p..</code> for integrating polynomials. For example the integral of <math>x^2</math> (which can be described in terms of its coefficients as <code>0 0 1</code>) is:
<lang j> 0 p.. 0 0 1
<syntaxhighlight lang="j"> 0 p.. 0 0 1
0 0 0 0.333333333333
0 0 0 0.333333333333
0 p.. 0 0 1x NB. or using rationals
0 p.. 0 0 1x NB. or using rationals
0 0 0 1r3</lang>
0 0 0 1r3</syntaxhighlight>
That is: <math>0x^0 + 0x^1 + 0x^2 + \tfrac{1}{3}x^3</math><br>
That is: <math>0x^0 + 0x^1 + 0x^2 + \tfrac{1}{3}x^3</math><br>
So to integrate <math>x^2</math> from 0 to &pi; :
So to integrate <math>x^2</math> from 0 to &pi; :
<lang j> 0 0 1 (0&p..@[ -~/@:p. ]) 0 1p1
<syntaxhighlight lang="j"> 0 0 1 (0&p..@[ -~/@:p. ]) 0 1p1
10.3354255601</lang>
10.3354255601</syntaxhighlight>


That said, J also has <code>d.</code> which can [http://www.jsoftware.com/help/dictionary/dddot.htm integrate] suitable functions.
That said, J also has <code>d.</code> which can [http://www.jsoftware.com/help/dictionary/dddot.htm integrate] suitable functions.


<lang j> *:d._1]1p1
<syntaxhighlight lang="j"> *:d._1]1p1
10.3354</lang>
10.3354</syntaxhighlight>


=={{header|Java}}==
=={{header|Java}}==
<lang java5>class NumericalIntegration
<syntaxhighlight lang="java5">class NumericalIntegration
{
{


Line 3,087: Line 3,087:
}
}
}
}
</syntaxhighlight>
</lang>


=={{header|Julia}}==
=={{header|Julia}}==
{{works with|Julia|0.6}}
{{works with|Julia|0.6}}


<lang julia>function simpson(f::Function, a::Number, b::Number, n::Integer)
<syntaxhighlight lang="julia">function simpson(f::Function, a::Number, b::Number, n::Integer)
h = (b - a) / n
h = (b - a) / n
s = f(a + h / 2)
s = f(a + h / 2)
Line 3,107: Line 3,107:
simpson(x -> x, 0, 6000, 6_000_000)
simpson(x -> x, 0, 6000, 6_000_000)


@show rst</lang>
@show rst</syntaxhighlight>


{{out}}
{{out}}
Line 3,113: Line 3,113:


=={{header|Kotlin}}==
=={{header|Kotlin}}==
<lang scala>// version 1.1.2
<syntaxhighlight lang="scala">// version 1.1.2


typealias Func = (Double) -> Double
typealias Func = (Double) -> Double
Line 3,138: Line 3,138:
integrate(0.0, 5000.0, 5_000_000) { it }
integrate(0.0, 5000.0, 5_000_000) { it }
integrate(0.0, 6000.0, 6_000_000) { it }
integrate(0.0, 6000.0, 6_000_000) { it }
}</lang>
}</syntaxhighlight>


{{out}}
{{out}}
Line 3,170: Line 3,170:
Following Python's presentation
Following Python's presentation


<syntaxhighlight lang="scheme">
<lang Scheme>
1) FUNCTIONS
1) FUNCTIONS


Line 3,258: Line 3,258:
trapezium 18006000
trapezium 18006000
simpson 18006000
simpson 18006000
</syntaxhighlight>
</lang>


=={{header|Liberty BASIC}}==
=={{header|Liberty BASIC}}==
Running the big loop value would take a VERY long time & seems unnecessary.<lang lb>
Running the big loop value would take a VERY long time & seems unnecessary.<syntaxhighlight lang="lb">
while 1
while 1
read x$
read x$
Line 3,379: Line 3,379:


end
end
</syntaxhighlight>
</lang>


Numerical integration
Numerical integration
Line 3,415: Line 3,415:


=={{header|Logo}}==
=={{header|Logo}}==
<lang logo>to i.left :fn :x :step
<syntaxhighlight lang="logo">to i.left :fn :x :step
output invoke :fn :x
output invoke :fn :x
end
end
Line 3,450: Line 3,450:
print integrate "i.mid "fn2 4 -1 2 ; 2.496091
print integrate "i.mid "fn2 4 -1 2 ; 2.496091
print integrate "i.trapezium "fn2 4 -1 2 ; 2.351014
print integrate "i.trapezium "fn2 4 -1 2 ; 2.351014
print integrate "i.simpsons "fn2 4 -1 2 ; 2.447732</lang>
print integrate "i.simpsons "fn2 4 -1 2 ; 2.447732</syntaxhighlight>


=={{header|Lua}}==
=={{header|Lua}}==
<lang lua>function leftRect( f, a, b, n )
<syntaxhighlight lang="lua">function leftRect( f, a, b, n )
local h = (b - a) / n
local h = (b - a) / n
local x = a
local x = a
Line 3,525: Line 3,525:
print( int_methods[i]( function(x) return x end, 0, 5000, 5000000 ) )
print( int_methods[i]( function(x) return x end, 0, 5000, 5000000 ) )
print( int_methods[i]( function(x) return x end, 0, 6000, 6000000 ) )
print( int_methods[i]( function(x) return x end, 0, 6000, 6000000 ) )
end</lang>
end</syntaxhighlight>


=={{header|Mathematica}}/{{header|Wolfram Language}}==
=={{header|Mathematica}}/{{header|Wolfram Language}}==
<lang Mathematica>leftRect[f_, a_Real, b_Real, N_Integer] :=
<syntaxhighlight lang="mathematica">leftRect[f_, a_Real, b_Real, N_Integer] :=
Module[{sum = 0, dx = (b - a)/N, x = a, n = N} ,
Module[{sum = 0, dx = (b - a)/N, x = a, n = N} ,
For[n = N, n > 0, n--, x += dx; sum += f[x];];
For[n = N, n > 0, n--, x += dx; sum += f[x];];
Line 3,553: Line 3,553:
For[n = 1, n < N, n++, sum1 += f[a + dx*n + dx/2];
For[n = 1, n < N, n++, sum1 += f[a + dx*n + dx/2];
sum2 += f[a + dx*n];];
sum2 += f[a + dx*n];];
Return [(dx/6)*(f[a] + f[b] + 4*sum1 + 2*sum2)]]</lang>
Return [(dx/6)*(f[a] + f[b] + 4*sum1 + 2*sum2)]]</syntaxhighlight>
<pre>f[x_] := x^3
<pre>f[x_] := x^3
g[x_] := 1/x
g[x_] := 1/x
Line 3,574: Line 3,574:


Function for performing left rectangular integration: leftRectIntegration.m
Function for performing left rectangular integration: leftRectIntegration.m
<lang MATLAB>function integral = leftRectIntegration(f,a,b,n)
<syntaxhighlight lang="matlab">function integral = leftRectIntegration(f,a,b,n)


format long;
format long;
Line 3,581: Line 3,581:
integral = width * sum( f(x(1:n-1)) );
integral = width * sum( f(x(1:n-1)) );
end</lang>
end</syntaxhighlight>


Function for performing right rectangular integration: rightRectIntegration.m
Function for performing right rectangular integration: rightRectIntegration.m
<lang MATLAB>function integral = rightRectIntegration(f,a,b,n)
<syntaxhighlight lang="matlab">function integral = rightRectIntegration(f,a,b,n)


format long;
format long;
Line 3,591: Line 3,591:
integral = width * sum( f(x(2:n)) );
integral = width * sum( f(x(2:n)) );
end</lang>
end</syntaxhighlight>


Function for performing mid-point rectangular integration: midPointRectIntegration.m
Function for performing mid-point rectangular integration: midPointRectIntegration.m
<lang MATLAB>function integral = midPointRectIntegration(f,a,b,n)
<syntaxhighlight lang="matlab">function integral = midPointRectIntegration(f,a,b,n)


format long;
format long;
Line 3,601: Line 3,601:
integral = width * sum( f( (x(1:n-1)+x(2:n))/2 ) );
integral = width * sum( f( (x(1:n-1)+x(2:n))/2 ) );
end</lang>
end</syntaxhighlight>


Function for performing trapezoidal integration: trapezoidalIntegration.m
Function for performing trapezoidal integration: trapezoidalIntegration.m
<lang MATLAB>function integral = trapezoidalIntegration(f,a,b,n)
<syntaxhighlight lang="matlab">function integral = trapezoidalIntegration(f,a,b,n)


format long;
format long;
Line 3,610: Line 3,610:
integral = trapz( x,f(x) );
integral = trapz( x,f(x) );
end</lang>
end</syntaxhighlight>


Simpson's rule for numerical integration is already included in Matlab as "quad()". It is not the same as the above examples, instead of specifying the amount of points to divide the x-axis into, the programmer passes the acceptable error tolerance for the calculation (parameter "tol").
Simpson's rule for numerical integration is already included in Matlab as "quad()". It is not the same as the above examples, instead of specifying the amount of points to divide the x-axis into, the programmer passes the acceptable error tolerance for the calculation (parameter "tol").
<lang MATLAB>integral = quad(f,a,b,tol)</lang>
<syntaxhighlight lang="matlab">integral = quad(f,a,b,tol)</syntaxhighlight>


Using anonymous functions
Using anonymous functions


<lang MATLAB>trapezoidalIntegration(@(x)( exp(-(x.^2)) ),0,10,100000)
<syntaxhighlight lang="matlab">trapezoidalIntegration(@(x)( exp(-(x.^2)) ),0,10,100000)


ans =
ans =


0.886226925452753</lang>
0.886226925452753</syntaxhighlight>


Using predefined functions
Using predefined functions


Built-in MATLAB function sin(x):
Built-in MATLAB function sin(x):
<lang MATLAB>quad(@sin,0,pi,1/1000000000000)
<syntaxhighlight lang="matlab">quad(@sin,0,pi,1/1000000000000)


ans =
ans =


2.000000000000000</lang>
2.000000000000000</syntaxhighlight>


User defined scripts and functions:
User defined scripts and functions:
fermiDirac.m
fermiDirac.m
<lang MATLAB>function answer = fermiDirac(x)
<syntaxhighlight lang="matlab">function answer = fermiDirac(x)
k = 8.617343e-5; %Boltazmann's Constant in eV/K
k = 8.617343e-5; %Boltazmann's Constant in eV/K
answer = 1./( 1+exp( (x)/(k*2000) ) ); %Fermi-Dirac distribution with mu = 0 and T = 2000K
answer = 1./( 1+exp( (x)/(k*2000) ) ); %Fermi-Dirac distribution with mu = 0 and T = 2000K
end</lang>
end</syntaxhighlight>


<lang MATLAB> rightRectIntegration(@fermiDirac,-1,1,1000000)
<syntaxhighlight lang="matlab"> rightRectIntegration(@fermiDirac,-1,1,1000000)


ans =
ans =


0.999998006023282</lang>
0.999998006023282</syntaxhighlight>


=={{header|Maxima}}==
=={{header|Maxima}}==
<lang maxima>right_rect(e, x, a, b, n) := block([h: (b - a) / n, s: 0],
<syntaxhighlight lang="maxima">right_rect(e, x, a, b, n) := block([h: (b - a) / n, s: 0],
for i from 1 thru n do s: s + subst(x = a + i * h, e),
for i from 1 thru n do s: s + subst(x = a + i * h, e),
s * h)$
s * h)$
Line 3,672: Line 3,672:
2 * log(2) - 1 - %, bfloat;
2 * log(2) - 1 - %, bfloat;


trapezium(1/x, x, 1, 100, 10000) - log(100), bfloat;</lang>
trapezium(1/x, x, 1, 100, 10000) - log(100), bfloat;</syntaxhighlight>


=={{header|Nim}}==
=={{header|Nim}}==
{{trans|Python}}
{{trans|Python}}
<lang nim>type Function = proc(x: float): float
<syntaxhighlight lang="nim">type Function = proc(x: float): float
type Rule = proc(f: Function; x, h: float): float
type Rule = proc(f: Function; x, h: float): float


Line 3,720: Line 3,720:
echo fName, " integrated using ", rName
echo fName, " integrated using ", rName
echo " from ", a, " to ", b, " (", steps, " steps) = ",
echo " from ", a, " to ", b, " (", steps, " steps) = ",
integrate(fun, float(a), float(b), steps, rule)</lang>
integrate(fun, float(a), float(b), steps, rule)</syntaxhighlight>


{{out}}
{{out}}
Line 3,766: Line 3,766:
=={{header|OCaml}}==
=={{header|OCaml}}==
The problem can be described as integrating using each of a set of methods, over a set of functions, so let us just build the solution in this modular way.
The problem can be described as integrating using each of a set of methods, over a set of functions, so let us just build the solution in this modular way.
First define the integration function:<lang ocaml>let integrate f a b steps meth =
First define the integration function:<syntaxhighlight lang="ocaml">let integrate f a b steps meth =
let h = (b -. a) /. float_of_int steps in
let h = (b -. a) /. float_of_int steps in
let rec helper i s =
let rec helper i s =
Line 3,772: Line 3,772:
else helper (succ i) (s +. meth f (a +. h *. float_of_int i) h)
else helper (succ i) (s +. meth f (a +. h *. float_of_int i) h)
in
in
h *. helper 0 0.</lang>Then list the methods:<lang ocaml>let methods = [
h *. helper 0 0.</syntaxhighlight>Then list the methods:<syntaxhighlight lang="ocaml">let methods = [
( "rect_l", fun f x _ -> f x);
( "rect_l", fun f x _ -> f x);
( "rect_m", fun f x h -> f (x +. h /. 2.) );
( "rect_m", fun f x h -> f (x +. h /. 2.) );
Line 3,778: Line 3,778:
( "trap", fun f x h -> (f x +. f (x +. h)) /. 2. );
( "trap", fun f x h -> (f x +. f (x +. h)) /. 2. );
( "simp", fun f x h -> (f x +. 4. *. f (x +. h /. 2.) +. f (x +. h)) /. 6. )
( "simp", fun f x h -> (f x +. 4. *. f (x +. h /. 2.) +. f (x +. h)) /. 6. )
]</lang> and functions (with limits and steps)<lang ocaml>let functions = [
]</syntaxhighlight> and functions (with limits and steps)<syntaxhighlight lang="ocaml">let functions = [
( "cubic", (fun x -> x*.x*.x), 0.0, 1.0, 100);
( "cubic", (fun x -> x*.x*.x), 0.0, 1.0, 100);
( "recip", (fun x -> 1.0/.x), 1.0, 100.0, 1000);
( "recip", (fun x -> 1.0/.x), 1.0, 100.0, 1000);
( "x to 5e3", (fun x -> x), 0.0, 5000.0, 5_000_000);
( "x to 5e3", (fun x -> x), 0.0, 5000.0, 5_000_000);
( "x to 6e3", (fun x -> x), 0.0, 6000.0, 6_000_000)
( "x to 6e3", (fun x -> x), 0.0, 6000.0, 6_000_000)
]</lang>and finally iterate the integration over both lists:<lang ocaml>let () =
]</syntaxhighlight>and finally iterate the integration over both lists:<syntaxhighlight lang="ocaml">let () =
List.iter (fun (s,f,lo,hi,n) ->
List.iter (fun (s,f,lo,hi,n) ->
Printf.printf "Testing function %s:\n" s;
Printf.printf "Testing function %s:\n" s;
Line 3,789: Line 3,789:
Printf.printf " method %s gives %.15g\n" name (integrate f lo hi n meth)
Printf.printf " method %s gives %.15g\n" name (integrate f lo hi n meth)
) methods
) methods
) functions</lang>Giving the output:
) functions</syntaxhighlight>Giving the output:
<pre>
<pre>
Testing function cubic:
Testing function cubic:
Line 3,819: Line 3,819:
=={{header|PARI/GP}}==
=={{header|PARI/GP}}==
Note also that double exponential integration is available as <code>intnum(x=a,b,f(x))</code> and Romberg integration is available as <code>intnumromb(x=a,b,f(x))</code>.
Note also that double exponential integration is available as <code>intnum(x=a,b,f(x))</code> and Romberg integration is available as <code>intnumromb(x=a,b,f(x))</code>.
<lang parigp>rectLeft(f, a, b, n)={
<syntaxhighlight lang="parigp">rectLeft(f, a, b, n)={
sum(i=0,n-1,f(a+(b-a)*i/n), 0.)*(b-a)/n
sum(i=0,n-1,f(a+(b-a)*i/n), 0.)*(b-a)/n
};
};
Line 3,847: Line 3,847:
test(x->1/x, 1, 100, 1000)
test(x->1/x, 1, 100, 1000)
test(x->x, 0, 5000, 5000000)
test(x->x, 0, 5000, 5000000)
test(x->x, 0, 6000, 6000000)</lang>
test(x->x, 0, 6000, 6000000)</syntaxhighlight>


Results:
Results:
Line 3,880: Line 3,880:


=={{header|Pascal}}==
=={{header|Pascal}}==
<lang pascal>function RectLeft(function f(x: real): real; xl, xr: real): real;
<syntaxhighlight lang="pascal">function RectLeft(function f(x: real): real; xl, xr: real): real;
begin
begin
RectLeft := f(xl)
RectLeft := f(xl)
Line 3,920: Line 3,920:
end;
end;
integrate := integral
integrate := integral
end;</lang>
end;</syntaxhighlight>


=={{header|Perl}}==
=={{header|Perl}}==
{{trans|Raku}}
{{trans|Raku}}
<lang perl>use feature 'say';
<syntaxhighlight lang="perl">use feature 'say';


sub leftrect {
sub leftrect {
Line 4,004: Line 4,004:
say for integrate('1 / $_', 1, 100, 1000, log(100)); say '';
say for integrate('1 / $_', 1, 100, 1000, log(100)); say '';
say for integrate('$_', 0, 5_000, 5_000_000, 12_500_000); say '';
say for integrate('$_', 0, 5_000, 5_000_000, 12_500_000); say '';
say for integrate('$_', 0, 6_000, 6_000_000, 18_000_000);</lang>
say for integrate('$_', 0, 6_000, 6_000_000, 18_000_000);</syntaxhighlight>
{{out}}
{{out}}
<pre>$_ ** 3
<pre>$_ ** 3
Line 4,043: Line 4,043:


=={{header|Phix}}==
=={{header|Phix}}==
<!--<lang Phix>(phixonline?)-->
<!--<syntaxhighlight lang="phix">(phixonline?)-->
<span style="color: #008080;">function</span> <span style="color: #000000;">rect_left</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">rid</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000080;font-style:italic;">/*h*/</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">rect_left</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">rid</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000080;font-style:italic;">/*h*/</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">rid</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">rid</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">)</span>
Line 4,116: Line 4,116:
<span style="color: #000000;">test</span><span style="color: #0000FF;">(</span><span style="color: #000000;">tests</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">test</span><span style="color: #0000FF;">(</span><span style="color: #000000;">tests</span><span style="color: #0000FF;">)</span>
<!--</lang>-->
<!--</syntaxhighlight>-->
{{out}}
{{out}}
<pre>
<pre>
Line 4,127: Line 4,127:


=={{header|PicoLisp}}==
=={{header|PicoLisp}}==
<lang PicoLisp>(scl 6)
<syntaxhighlight lang="picolisp">(scl 6)


(de leftRect (Fun X)
(de leftRect (Fun X)
Line 4,158: Line 4,158:
(*/ H Sum 1.0) ) )
(*/ H Sum 1.0) ) )


(prinl (round (integrate square 3.0 7.0 30 simpson)))</lang>
(prinl (round (integrate square 3.0 7.0 30 simpson)))</syntaxhighlight>
Output:
Output:
<pre>105.333</pre>
<pre>105.333</pre>


=={{header|PL/I}}==
=={{header|PL/I}}==
<lang PL/I>integrals: procedure options (main); /* 1 September 2019 */
<syntaxhighlight lang="pl/i">integrals: procedure options (main); /* 1 September 2019 */


f: procedure (x, function) returns (float(18));
f: procedure (x, function) returns (float(18));
Line 4,234: Line 4,234:


end integrals;
end integrals;
</syntaxhighlight>
</lang>
<pre>
<pre>
Rectangle-left Rectangle-mid Rectangle-right Trapezoid Simpson
Rectangle-left Rectangle-mid Rectangle-right Trapezoid Simpson
Line 4,245: Line 4,245:
=={{header|PureBasic}}==
=={{header|PureBasic}}==


<lang PureBasic>Prototype.d TestFunction(Arg.d)
<syntaxhighlight lang="purebasic">Prototype.d TestFunction(Arg.d)


Procedure.d LeftIntegral(Start, Stop, Steps, *func.TestFunction)
Procedure.d LeftIntegral(Start, Stop, Steps, *func.TestFunction)
Line 4,346: Line 4,346:
Answer$+"Trapezium="+StrD(Trapezium (0,6000,6000000,@Test3()))+#CRLF$
Answer$+"Trapezium="+StrD(Trapezium (0,6000,6000000,@Test3()))+#CRLF$
Answer$+"Simpson ="+StrD(Simpson (0,6000,6000000,@Test3()))
Answer$+"Simpson ="+StrD(Simpson (0,6000,6000000,@Test3()))
MessageRequester("Answer should be 18,000,000",Answer$) </lang>
MessageRequester("Answer should be 18,000,000",Answer$) </syntaxhighlight>
<pre>Left =0.2353220100
<pre>Left =0.2353220100
Mid =0.2401367513
Mid =0.2401367513
Line 4,373: Line 4,373:
=={{header|Python}}==
=={{header|Python}}==
Answers are first given using floating point arithmatic, then using fractions, only converted to floating point on output.
Answers are first given using floating point arithmatic, then using fractions, only converted to floating point on output.
<lang python>from fractions import Fraction
<syntaxhighlight lang="python">from fractions import Fraction


def left_rect(f,x,h):
def left_rect(f,x,h):
Line 4,427: Line 4,427:
print('%s integrated using %s\n from %r to %r (%i steps and fractions) = %r' %
print('%s integrated using %s\n from %r to %r (%i steps and fractions) = %r' %
(func.__name__, rule.__name__, a, b, steps,
(func.__name__, rule.__name__, a, b, steps,
float(integrate( func, a, b, steps, rule))))</lang>
float(integrate( func, a, b, steps, rule))))</syntaxhighlight>


'''Tests'''
'''Tests'''
<lang python>for a, b, steps, func in ((0., 1., 100, cube), (1., 100., 1000, reciprocal)):
<syntaxhighlight lang="python">for a, b, steps, func in ((0., 1., 100, cube), (1., 100., 1000, reciprocal)):
for rule in (left_rect, mid_rect, right_rect, trapezium, simpson):
for rule in (left_rect, mid_rect, right_rect, trapezium, simpson):
print('%s integrated using %s\n from %r to %r (%i steps) = %r' %
print('%s integrated using %s\n from %r to %r (%i steps) = %r' %
Line 4,452: Line 4,452:
print('%s integrated using %s\n from %r to %r (%i steps and fractions) = %r' %
print('%s integrated using %s\n from %r to %r (%i steps and fractions) = %r' %
(func.__name__, rule.__name__, a, b, steps,
(func.__name__, rule.__name__, a, b, steps,
float(integrate( func, a, b, steps, rule))))</lang>
float(integrate( func, a, b, steps, rule))))</syntaxhighlight>


'''Sample test Output'''
'''Sample test Output'''
Line 4,537: Line 4,537:


A faster Simpson's rule integrator is
A faster Simpson's rule integrator is
<lang python>def faster_simpson(f, a, b, steps):
<syntaxhighlight lang="python">def faster_simpson(f, a, b, steps):
h = (b-a)/float(steps)
h = (b-a)/float(steps)
a1 = a+h/2
a1 = a+h/2
s1 = sum( f(a1+i*h) for i in range(0,steps))
s1 = sum( f(a1+i*h) for i in range(0,steps))
s2 = sum( f(a+i*h) for i in range(1,steps))
s2 = sum( f(a+i*h) for i in range(1,steps))
return (h/6.0)*(f(a)+f(b)+4.0*s1+2.0*s2)</lang>
return (h/6.0)*(f(a)+f(b)+4.0*s1+2.0*s2)</syntaxhighlight>


=={{header|R}}==
=={{header|R}}==
The integ function defined below uses arbitrary abscissae and weights passed as argument (resp. u and v). It assumes that f can take a vector argument.
The integ function defined below uses arbitrary abscissae and weights passed as argument (resp. u and v). It assumes that f can take a vector argument.


<lang rsplus>integ <- function(f, a, b, n, u, v) {
<syntaxhighlight lang="rsplus">integ <- function(f, a, b, n, u, v) {
h <- (b - a) / n
h <- (b - a) / n
s <- 0
s <- 0
Line 4,578: Line 4,578:
test(\(x) x, 0, 6000, 6e6)
test(\(x) x, 0, 6000, 6e6)
# rect.left rect.right rect.mid trapezoidal simpson
# rect.left rect.right rect.mid trapezoidal simpson
# 1.8e+07 1.8e+07 1.8e+07 1.8e+07 1.8e+07</lang>
# 1.8e+07 1.8e+07 1.8e+07 1.8e+07 1.8e+07</syntaxhighlight>


=={{header|Racket}}==
=={{header|Racket}}==
<lang racket>
<syntaxhighlight lang="racket">
#lang racket
#lang racket
(define (integrate f a b steps meth)
(define (integrate f a b steps meth)
Line 4,605: Line 4,605:
(test (λ(x) x) 0. 5000. 5000000 "IDENTITY")
(test (λ(x) x) 0. 5000. 5000000 "IDENTITY")
(test (λ(x) x) 0. 6000. 6000000 "IDENTITY")
(test (λ(x) x) 0. 6000. 6000000 "IDENTITY")
</syntaxhighlight>
</lang>
Output:
Output:
<lang racket>
<syntaxhighlight lang="racket">
CUBED
CUBED
left-rect: 0.24502500000000005
left-rect: 0.24502500000000005
Line 4,635: Line 4,635:
trapezium: 17999999.999999993
trapezium: 17999999.999999993
simpson: 17999999.999999993
simpson: 17999999.999999993
</syntaxhighlight>
</lang>


=={{header|Raku}}==
=={{header|Raku}}==
Line 4,644: Line 4,644:
{{works with|Rakudo|2018.09}}
{{works with|Rakudo|2018.09}}


<lang perl6>use MONKEY-SEE-NO-EVAL;
<syntaxhighlight lang="raku" line>use MONKEY-SEE-NO-EVAL;


sub leftrect(&f, $a, $b, $n) {
sub leftrect(&f, $a, $b, $n) {
Line 4,711: Line 4,711:
.say for integrate '1 / *', 1, 100, 1000, log(100); say '';
.say for integrate '1 / *', 1, 100, 1000, log(100); say '';
.say for integrate '*.self', 0, 5_000, 5_000_000, 12_500_000; say '';
.say for integrate '*.self', 0, 5_000, 5_000_000, 12_500_000; say '';
.say for integrate '*.self', 0, 6_000, 6_000_000, 18_000_000;</lang>
.say for integrate '*.self', 0, 6_000, 6_000_000, 18_000_000;</syntaxhighlight>
{{out}}
{{out}}
<pre>{ $_ ** 3 }
<pre>{ $_ ** 3 }
Line 4,751: Line 4,751:
=={{header|REXX}}==
=={{header|REXX}}==
Note: &nbsp; there was virtually no difference in accuracy between &nbsp; '''numeric digits 9''' &nbsp; (the default) &nbsp; and &nbsp; '''numeric digits 20'''.
Note: &nbsp; there was virtually no difference in accuracy between &nbsp; '''numeric digits 9''' &nbsp; (the default) &nbsp; and &nbsp; '''numeric digits 20'''.
<lang rexx>/*REXX pgm performs numerical integration using 5 different algorithms and show results.*/
<syntaxhighlight lang="rexx">/*REXX pgm performs numerical integration using 5 different algorithms and show results.*/
numeric digits 20 /*use twenty decimal digits precision. */
numeric digits 20 /*use twenty decimal digits precision. */


Line 4,798: Line 4,798:
do x=a by h for #; $= $ + (f(x) + f(x+h))
do x=a by h for #; $= $ + (f(x) + f(x+h))
end /*x*/
end /*x*/
return $*h/2</lang>
return $*h/2</syntaxhighlight>
{{out|output|text=&nbsp; when using the default inputs:}}
{{out|output|text=&nbsp; when using the default inputs:}}
<pre>
<pre>
Line 4,831: Line 4,831:


=={{header|Ring}}==
=={{header|Ring}}==
<lang ring>
<syntaxhighlight lang="ring">
# Project : Numerical integration
# Project : Numerical integration


Line 4,917: Line 4,917:
eval("result = " + x2)
eval("result = " + x2)
return (d / 6) * (f + result + 4 * s1 + 2 * s)
return (d / 6) * (f + result + 4 * s1 + 2 * s)
</syntaxhighlight>
</lang>
Output:
Output:
<pre>
<pre>
Line 4,929: Line 4,929:
=={{header|Ruby}}==
=={{header|Ruby}}==
{{trans|Tcl}}
{{trans|Tcl}}
<lang ruby>def leftrect(f, left, right)
<syntaxhighlight lang="ruby">def leftrect(f, left, right)
f.call(left)
f.call(left)
end
end
Line 4,986: Line 4,986:
printf " %-10s %s\t(%.1f%%)\n", method, int, diff
printf " %-10s %s\t(%.1f%%)\n", method, int, diff
end
end
end</lang>
end</syntaxhighlight>
outputs
outputs
<pre>integral of #<Method: Object#square> from 0 to 3.14159265358979 in 10 steps
<pre>integral of #<Method: Object#square> from 0 to 3.14159265358979 in 10 steps
Line 5,003: Line 5,003:
=={{header|Rust}}==
=={{header|Rust}}==
This is a partial solution and only implements trapezium integration.
This is a partial solution and only implements trapezium integration.
<lang rust>fn integral<F>(f: F, range: std::ops::Range<f64>, n_steps: u32) -> f64
<syntaxhighlight lang="rust">fn integral<F>(f: F, range: std::ops::Range<f64>, n_steps: u32) -> f64
where F: Fn(f64) -> f64
where F: Fn(f64) -> f64
{
{
Line 5,022: Line 5,022:
println!("{}", integral(|x| x, 0.0..5000.0, 5_000_000));
println!("{}", integral(|x| x, 0.0..5000.0, 5_000_000));
println!("{}", integral(|x| x, 0.0..6000.0, 6_000_000));
println!("{}", integral(|x| x, 0.0..6000.0, 6_000_000));
}</lang>
}</syntaxhighlight>


{{out}}
{{out}}
Line 5,031: Line 5,031:


=={{header|Scala}}==
=={{header|Scala}}==
<lang scala>object NumericalIntegration {
<syntaxhighlight lang="scala">object NumericalIntegration {
def leftRect(f:Double=>Double, a:Double, b:Double)=f(a)
def leftRect(f:Double=>Double, a:Double, b:Double)=f(a)
def midRect(f:Double=>Double, a:Double, b:Double)=f((a+b)/2)
def midRect(f:Double=>Double, a:Double, b:Double)=f((a+b)/2)
Line 5,065: Line 5,065:
print(fn3, 0, 6000, 6000000)
print(fn3, 0, 6000, 6000000)
}
}
}</lang>
}</syntaxhighlight>
Output:
Output:
<pre>rectangular left : 0,245025
<pre>rectangular left : 0,245025
Line 5,093: Line 5,093:
=={{header|Scheme}}==
=={{header|Scheme}}==


<lang scheme>(define (integrate f a b steps meth)
<syntaxhighlight lang="scheme">(define (integrate f a b steps meth)
(define h (/ (- b a) steps))
(define h (/ (- b a) steps))
(* h
(* h
Line 5,113: Line 5,113:
(define rr (integrate square 0 1 10 right-rect))
(define rr (integrate square 0 1 10 right-rect))
(define t (integrate square 0 1 10 trapezium))
(define t (integrate square 0 1 10 trapezium))
(define s (integrate square 0 1 10 simpson))</lang>
(define s (integrate square 0 1 10 simpson))</syntaxhighlight>


=={{header|SequenceL}}==
=={{header|SequenceL}}==
<lang sequencel>import <Utilities/Conversion.sl>;
<syntaxhighlight lang="sequencel">import <Utilities/Conversion.sl>;
import <Utilities/Sequence.sl>;
import <Utilities/Sequence.sl>;


Line 5,175: Line 5,175:
delimit(delimit(heading ++ transpose(funcs ++ ranges ++ trimEndZeroes(floatToString(tests, 8))), '\t'), '\n');
delimit(delimit(heading ++ transpose(funcs ++ ranges ++ trimEndZeroes(floatToString(tests, 8))), '\t'), '\n');


trimEndZeroes(x(1)) := x when size(x) = 0 else x when x[size(x)] /= '0' else trimEndZeroes(x[1...size(x)-1]);</lang>
trimEndZeroes(x(1)) := x when size(x) = 0 else x when x[size(x)] /= '0' else trimEndZeroes(x[1...size(x)-1]);</syntaxhighlight>


{{out}}
{{out}}
Line 5,188: Line 5,188:
=={{header|Sidef}}==
=={{header|Sidef}}==
{{trans|Raku}}
{{trans|Raku}}
<lang ruby>func sum(f, start, from, to) {
<syntaxhighlight lang="ruby">func sum(f, start, from, to) {
var s = 0;
var s = 0;
RangeNum(start, to, from-start).each { |i|
RangeNum(start, to, from-start).each { |i|
Line 5,241: Line 5,241:
tryem('1/x', { 1 / _ }, 1, 100, 1000, log(100));
tryem('1/x', { 1 / _ }, 1, 100, 1000, log(100));
tryem('x', { _ }, 0, 5_000, 5_000_000, 12_500_000);
tryem('x', { _ }, 0, 5_000, 5_000_000, 12_500_000);
tryem('x', { _ }, 0, 6_000, 6_000_000, 18_000_000);</lang>
tryem('x', { _ }, 0, 6_000, 6_000_000, 18_000_000);</syntaxhighlight>


=={{header|Standard ML}}==
=={{header|Standard ML}}==
<lang sml>fun integrate (f, a, b, steps, meth) = let
<syntaxhighlight lang="sml">fun integrate (f, a, b, steps, meth) = let
val h = (b - a) / real steps
val h = (b - a) / real steps
fun helper (i, s) =
fun helper (i, s) =
Line 5,266: Line 5,266:
val rr = integrate (square, 0.0, 1.0, 10, right_rect)
val rr = integrate (square, 0.0, 1.0, 10, right_rect)
val t = integrate (square, 0.0, 1.0, 10, trapezium )
val t = integrate (square, 0.0, 1.0, 10, trapezium )
val s = integrate (square, 0.0, 1.0, 10, simpson )</lang>
val s = integrate (square, 0.0, 1.0, 10, simpson )</syntaxhighlight>


=={{header|Stata}}==
=={{header|Stata}}==
<lang>mata
<syntaxhighlight lang="text">mata
function integrate(f,a,b,n,u,v) {
function integrate(f,a,b,n,u,v) {
s = 0
s = 0
Line 5,309: Line 5,309:
test(&id(),0,5000,5000000)
test(&id(),0,5000,5000000)
test(&id(),0,6000,6000000)
test(&id(),0,6000,6000000)
end</lang>
end</syntaxhighlight>


'''Output'''
'''Output'''
Line 5,334: Line 5,334:


=={{header|Swift}}==
=={{header|Swift}}==
<lang swift>public enum IntegrationType : CaseIterable {
<syntaxhighlight lang="swift">public enum IntegrationType : CaseIterable {
case rectangularLeft
case rectangularLeft
case rectangularRight
case rectangularRight
Line 5,438: Line 5,438:
print("f(x) = 1 / x:", types.map({ integrate(from: 1, to: 100, n: 1000, using: $0, f: { 1 / $0 }) }))
print("f(x) = 1 / x:", types.map({ integrate(from: 1, to: 100, n: 1000, using: $0, f: { 1 / $0 }) }))
print("f(x) = x, 0 -> 5_000:", types.map({ integrate(from: 0, to: 5_000, n: 5_000_000, using: $0, f: { $0 }) }))
print("f(x) = x, 0 -> 5_000:", types.map({ integrate(from: 0, to: 5_000, n: 5_000_000, using: $0, f: { $0 }) }))
print("f(x) = x, 0 -> 6_000:", types.map({ integrate(from: 0, to: 6_000, n: 6_000_000, using: $0, f: { $0 }) }))</lang>
print("f(x) = x, 0 -> 6_000:", types.map({ integrate(from: 0, to: 6_000, n: 6_000_000, using: $0, f: { $0 }) }))</syntaxhighlight>


{{out}}
{{out}}
Line 5,447: Line 5,447:


=={{header|Tcl}}==
=={{header|Tcl}}==
<lang tcl>package require Tcl 8.5
<syntaxhighlight lang="tcl">package require Tcl 8.5


proc leftrect {f left right} {
proc leftrect {f left right} {
Line 5,500: Line 5,500:
puts [format " %-10s %s\t(%.1f%%)" $method $int $diff]
puts [format " %-10s %s\t(%.1f%%)" $method $int $diff]
}
}
}</lang>
}</syntaxhighlight>
<pre>integral of square(x) from 0 to 3.141592653589793 in 10 steps
<pre>integral of square(x) from 0 to 3.141592653589793 in 10 steps
leftrect 8.836788853885448 (-14.5%)
leftrect 8.836788853885448 (-14.5%)
Line 5,530: Line 5,530:
the integrand <math>f</math>, the bounds <math>(a,b)</math>, and the number of intervals <math>n</math>.
the integrand <math>f</math>, the bounds <math>(a,b)</math>, and the number of intervals <math>n</math>.


<lang Ursala>#import std
<syntaxhighlight lang="ursala">#import std
#import nat
#import nat
#import flo
#import flo
Line 5,536: Line 5,536:
(integral_by "m") ("f","a","b","n") =
(integral_by "m") ("f","a","b","n") =


iprod ^(* ! div\float"n" minus/"b" "a",~&) ("m" "f")*ytp (ari successor "n")/"a" "b"</lang>
iprod ^(* ! div\float"n" minus/"b" "a",~&) ("m" "f")*ytp (ari successor "n")/"a" "b"</syntaxhighlight>
An alternative way of defining this function shown below prevents redundant evaluations of the integrand
An alternative way of defining this function shown below prevents redundant evaluations of the integrand
at the cost of building a table-driven finite map in advance.
at the cost of building a table-driven finite map in advance.
<lang Ursala>(integral_by "m") ("f","a","b","n") =
<syntaxhighlight lang="ursala">(integral_by "m") ("f","a","b","n") =


iprod ^(* ! div\float"n" minus/"b" "a",~&) ^H(*+ "m"+ -:"f"+ * ^/~& "f",~&ytp) (ari successor "n")/"a" "b"</lang>
iprod ^(* ! div\float"n" minus/"b" "a",~&) ^H(*+ "m"+ -:"f"+ * ^/~& "f",~&ytp) (ari successor "n")/"a" "b"</syntaxhighlight>
As mentioned in the Haskell solution, the latter choice is preferable if evaluating the integrand
As mentioned in the Haskell solution, the latter choice is preferable if evaluating the integrand
is expensive.
is expensive.
An integrating function is defined for each method as follows.
An integrating function is defined for each method as follows.
<lang Ursala>left = integral_by "f". ("l","r"). "f" "l"
<syntaxhighlight lang="ursala">left = integral_by "f". ("l","r"). "f" "l"
right = integral_by "f". ("l","r"). "f" "r"
right = integral_by "f". ("l","r"). "f" "r"
midpoint = integral_by "f". ("l","r"). "f" div\2. plus/"l" "r"
midpoint = integral_by "f". ("l","r"). "f" div\2. plus/"l" "r"
trapezium = integral_by "f". ("l","r"). div\2. plus "f"~~/"l" "r"
trapezium = integral_by "f". ("l","r"). div\2. plus "f"~~/"l" "r"
simpson = integral_by "f". ("l","r"). div\6. plus:-0. <"f" "l",times/4. "f" div\2. plus/"l" "r","f" "r"></lang>
simpson = integral_by "f". ("l","r"). div\6. plus:-0. <"f" "l",times/4. "f" div\2. plus/"l" "r","f" "r"></syntaxhighlight>
As shown above, the method passed to the <code>integral_by</code> function
As shown above, the method passed to the <code>integral_by</code> function
is itself a higher order function taking an integrand <math>f</math> as an argument and
is itself a higher order function taking an integrand <math>f</math> as an argument and
Line 5,555: Line 5,555:
Here is a test program showing the results of integrating the square from zero to <math>\pi</math> in ten intervals
Here is a test program showing the results of integrating the square from zero to <math>\pi</math> in ten intervals
by all five methods.
by all five methods.
<lang Ursala>#cast %eL
<syntaxhighlight lang="ursala">#cast %eL


examples = <.left,midpoint,rignt,trapezium,simpson> (sqr,0.,pi,10)</lang>
examples = <.left,midpoint,rignt,trapezium,simpson> (sqr,0.,pi,10)</syntaxhighlight>
output:
output:
<pre>
<pre>
Line 5,573: Line 5,573:
The following program does not follow the task requirement on two points: first, the same function is used for all quadrature methods, as they are really the same thing with different parameters (abscissas and weights). And since it's getting rather slow for a large number of intervals, the last two are integrated with resp. 50,000 and 60,000 intervals. It does not make sense anyway to use more, for such a simple function (and if really it were difficult to integrate, one would rely one more sophistcated methods).
The following program does not follow the task requirement on two points: first, the same function is used for all quadrature methods, as they are really the same thing with different parameters (abscissas and weights). And since it's getting rather slow for a large number of intervals, the last two are integrated with resp. 50,000 and 60,000 intervals. It does not make sense anyway to use more, for such a simple function (and if really it were difficult to integrate, one would rely one more sophistcated methods).


<lang vb>Option Explicit
<syntaxhighlight lang="vb">Option Explicit
Option Base 1
Option Base 1


Line 5,628: Line 5,628:
Next j
Next j
Next i
Next i
End Sub</lang>
End Sub</syntaxhighlight>


=={{header|Wren}}==
=={{header|Wren}}==
{{trans|Kotlin}}
{{trans|Kotlin}}
{{libheader|Wren-fmt}}
{{libheader|Wren-fmt}}
<lang ecmascript>import "/fmt" for Fmt
<syntaxhighlight lang="ecmascript">import "/fmt" for Fmt


var integrate = Fn.new { |a, b, n, f|
var integrate = Fn.new { |a, b, n, f|
Line 5,655: Line 5,655:
integrate.call(0, 5000, 5000000) { |v| v }
integrate.call(0, 5000, 5000000) { |v| v }
integrate.call(0, 6000, 6000000) { |v| v }
integrate.call(0, 6000, 6000000) { |v| v }
</syntaxhighlight>
</lang>


{{out}}
{{out}}
Line 5,685: Line 5,685:


=={{header|XPL0}}==
=={{header|XPL0}}==
<lang XPL0>include c:\cxpl\codes; \intrinsic 'code' declarations
<syntaxhighlight lang="xpl0">include c:\cxpl\codes; \intrinsic 'code' declarations


func real Func(FN, X); \Return F(X) for function number FN
func real Func(FN, X); \Return F(X) for function number FN
Line 5,731: Line 5,731:
Integrate(0.0, 5000.0, 3, 5_000_000);
Integrate(0.0, 5000.0, 3, 5_000_000);
Integrate(0.0, 6000.0, 3, 6_000_000);
Integrate(0.0, 6000.0, 3, 6_000_000);
]</lang>
]</syntaxhighlight>


Interestingly, the small rounding errors creep in when millions of
Interestingly, the small rounding errors creep in when millions of
Line 5,747: Line 5,747:
=={{header|Yabasic}}==
=={{header|Yabasic}}==
Based on the XPL0entry and the Free BASIC entry
Based on the XPL0entry and the Free BASIC entry
<lang yabasic>// Rosetta Code problem: https://rosettacode.org/wiki/Numerical_integration
<syntaxhighlight lang="yabasic">// Rosetta Code problem: https://rosettacode.org/wiki/Numerical_integration
// by Jjuanhdez, 06/2022
// by Jjuanhdez, 06/2022


Line 5,819: Line 5,819:
print " ";
print " ";
print str$(Area, frmt$)
print str$(Area, frmt$)
end sub</lang>
end sub</syntaxhighlight>


=={{header|zkl}}==
=={{header|zkl}}==
{{trans|D}}
{{trans|D}}
<lang zkl>fcn integrate(F,f,a,b,steps){
<syntaxhighlight lang="zkl">fcn integrate(F,f,a,b,steps){
h:=(b - a) / steps;
h:=(b - a) / steps;
h*(0).reduce(steps,'wrap(s,i){ F(f, h*i + a, h) + s },0.0);
h*(0).reduce(steps,'wrap(s,i){ F(f, h*i + a, h) + s },0.0);
Line 5,846: Line 5,846:
"%s %f".fmt(nm,integrate(f,a.xplode())).println() }, fs);
"%s %f".fmt(nm,integrate(f,a.xplode())).println() }, fs);
println();
println();
}</lang>
}</syntaxhighlight>
{{out}}
{{out}}
<pre>
<pre>