Numerical integration: Difference between revisions
Content added Content deleted
m (→{{header|R}}) |
Thundergnat (talk | contribs) m (syntax highlighting fixup automation) |
||
Line 48: | Line 48: | ||
{{trans|Nim}} |
{{trans|Nim}} |
||
< |
<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)))</ |
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: |
||
< |
<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); |
||
}</ |
}</syntaxhighlight> |
||
Usage: |
Usage: |
||
< |
<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: |
||
< |
<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;</ |
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: |
||
< |
<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;</ |
end Integrate;</syntaxhighlight> |
||
Test driver: |
Test driver: |
||
< |
<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}}== |
||
< |
<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</ |
SKIP</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 513: | Line 513: | ||
=={{header|ALGOL W}}== |
=={{header|ALGOL W}}== |
||
{{Trans|ALGOL 68}} |
{{Trans|ALGOL 68}} |
||
< |
<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.</ |
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] |
||
< |
<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 |
||
}</ |
}</syntaxhighlight> |
||
=={{header|BASIC}}== |
=={{header|BASIC}}== |
||
{{works with|QuickBasic|4.5}} |
{{works with|QuickBasic|4.5}} |
||
{{trans|Java}} |
{{trans|Java}} |
||
< |
<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</ |
END FUNCTION</syntaxhighlight> |
||
=={{header|BBC BASIC}}== |
=={{header|BBC BASIC}}== |
||
< |
<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)</ |
= (d / 6) * (f + EVAL(x$) + 4 * s1 + 2 * s2)</syntaxhighlight> |
||
'''Output:''' |
'''Output:''' |
||
<pre> |
<pre> |
||
Line 818: | Line 818: | ||
=={{header|C}}== |
=={{header|C}}== |
||
< |
<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); |
||
}</ |
}</syntaxhighlight> |
||
< |
<syntaxhighlight lang="c">/* test */ |
||
double f3(double x) |
double f3(double x) |
||
{ |
{ |
||
Line 947: | Line 947: | ||
printf("\n"); |
printf("\n"); |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
=={{header|C sharp|C#}}== |
=={{header|C sharp|C#}}== |
||
< |
<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); |
||
} |
} |
||
}</ |
}</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</ |
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. |
||
< |
<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());</ |
double s = integrate(f, 0.0, 1.0, 10, simpson());</syntaxhighlight> |
||
=={{header|Chapel}}== |
=={{header|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}} |
||
< |
<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}}== |
||
< |
<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))))))</ |
(* 2 sum2))))))</syntaxhighlight> |
||
=={{header|D}}== |
=={{header|D}}== |
||
< |
<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(); |
||
} |
} |
||
}</ |
}</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: |
||
< |
<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(); |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
=={{header|Delphi}}== |
=={{header|Delphi}}== |
||
{{libheader| System.SysUtils}} |
{{libheader| System.SysUtils}} |
||
{{Trans|Python}} |
{{Trans|Python}} |
||
< |
<syntaxhighlight lang="delphi">program Numerical_integration; |
||
{$APPTYPE CONSOLE} |
{$APPTYPE CONSOLE} |
||
Line 1,785: | Line 1,785: | ||
end; |
end; |
||
readln; |
readln; |
||
end.</ |
end.</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>Integrate x^3 |
<pre>Integrate x^3 |
||
Line 1,815: | Line 1,815: | ||
{{trans|Python}} |
{{trans|Python}} |
||
< |
<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) } |
||
}</ |
}</syntaxhighlight> |
||
< |
<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</ |
# value: 0.10000000002160479</syntaxhighlight> |
||
=={{header|Elixir}}== |
=={{header|Elixir}}== |
||
< |
<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)</ |
Numerical.integrate(f4, 0, 6000, 6_000_000)</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 1,918: | Line 1,918: | ||
=={{header|Euphoria}}== |
=={{header|Euphoria}}== |
||
< |
<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"))</ |
? int_simpson({0,10},1000,routine_id("x"))</syntaxhighlight> |
||
Output: |
Output: |
||
Line 2,016: | Line 2,016: | ||
=={{header|F Sharp}}== |
=={{header|F Sharp}}== |
||
< |
<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}}== |
||
< |
<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</ |
18000000</syntaxhighlight> |
||
=={{header|Forth}}== |
=={{header|Forth}}== |
||
< |
<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</ |
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: |
||
< |
<syntaxhighlight lang="fortran">elemental function elemf(x) |
||
real :: elemf, x |
real :: elemf, x |
||
elemf = f(x) |
elemf = f(x) |
||
end function elemf</ |
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. |
||
< |
<syntaxhighlight lang="fortran">module Integration |
||
implicit none |
implicit none |
||
Line 2,183: | Line 2,183: | ||
end function integrate |
end function integrate |
||
end module Integration</ |
end module Integration</syntaxhighlight> |
||
Usage example: |
Usage example: |
||
< |
<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</ |
end program IntegrationTest</syntaxhighlight> |
||
The FunctionHolder module: |
The FunctionHolder module: |
||
< |
<syntaxhighlight lang="fortran">module FunctionHolder |
||
implicit none |
implicit none |
||
Line 2,213: | Line 2,213: | ||
end function afun |
end function afun |
||
end module FunctionHolder</ |
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 |
||
< |
<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</ |
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}}== |
||
< |
<syntaxhighlight lang="go">package main |
||
import ( |
import ( |
||
Line 2,513: | Line 2,513: | ||
fmt.Println("") |
fmt.Println("") |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 2,554: | Line 2,554: | ||
=={{header|Groovy}}== |
=={{header|Groovy}}== |
||
Solution: |
Solution: |
||
< |
<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())) |
||
} |
} |
||
}</ |
}</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). |
||
< |
<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</ |
assert ((simpsonsIntegral([-Math.E, Math.PI], 1, cubicPoly) - cpIntegralCalcMinusEToPI)/ cpIntegralCalcMinusEToPI).abs() < tolerance**2.5 // 1 in 10 billion</syntaxhighlight> |
||
Requested Demonstrations: |
Requested Demonstrations: |
||
< |
<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 ()</ |
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 |
||
< |
<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: |
||
< |
<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]]</ |
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: |
||
< |
<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</ |
inter n (y:ys) = y : inter n ys</syntaxhighlight> |
||
And now we can just define |
And now we can just define |
||
< |
<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]</ |
intSimpson = integrateClosed 3 [1,4,1]</syntaxhighlight> |
||
or, as easily, some additional schemes: |
or, as easily, some additional schemes: |
||
< |
<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]</ |
intOpen2 = integrateOpen 3 [8,-4,8]</syntaxhighlight> |
||
Some examples: |
Some examples: |
||
Line 2,769: | Line 2,769: | ||
The whole program: |
The whole program: |
||
< |
<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 ' ')</ |
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:=== |
||
< |
<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'</ |
simpson =: adverb def '6 %~ +/ 1 1 4 * u y, -:+/y'</syntaxhighlight> |
||
===Example usage=== |
===Example usage=== |
||
====Required Examples==== |
====Required Examples==== |
||
< |
<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</ |
1.8e7</syntaxhighlight> |
||
====Older Examples==== |
====Older Examples==== |
||
Integrate <code>square</code> (<code>*:</code>) from 0 to π in 10 steps using various methods. |
Integrate <code>square</code> (<code>*:</code>) from 0 to π in 10 steps using various methods. |
||
< |
<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</ |
10.3354255601</syntaxhighlight> |
||
Integrate <code>sin</code> from 0 to π in 10 steps using various methods. |
Integrate <code>sin</code> from 0 to π in 10 steps using various methods. |
||
< |
<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</ |
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: |
||
< |
<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</ |
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 π : |
So to integrate <math>x^2</math> from 0 to π : |
||
< |
<syntaxhighlight lang="j"> 0 0 1 (0&p..@[ -~/@:p. ]) 0 1p1 |
||
10.3354255601</ |
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. |
||
< |
<syntaxhighlight lang="j"> *:d._1]1p1 |
||
10.3354</ |
10.3354</syntaxhighlight> |
||
=={{header|Java}}== |
=={{header|Java}}== |
||
< |
<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}} |
||
< |
<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</ |
@show rst</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 3,113: | Line 3,113: | ||
=={{header|Kotlin}}== |
=={{header|Kotlin}}== |
||
< |
<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 } |
||
}</ |
}</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.< |
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}}== |
||
< |
<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</ |
print integrate "i.simpsons "fn2 4 -1 2 ; 2.447732</syntaxhighlight> |
||
=={{header|Lua}}== |
=={{header|Lua}}== |
||
< |
<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</ |
end</syntaxhighlight> |
||
=={{header|Mathematica}}/{{header|Wolfram Language}}== |
=={{header|Mathematica}}/{{header|Wolfram Language}}== |
||
< |
<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)]]</ |
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 |
||
< |
<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</ |
end</syntaxhighlight> |
||
Function for performing right rectangular integration: rightRectIntegration.m |
Function for performing right rectangular integration: rightRectIntegration.m |
||
< |
<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</ |
end</syntaxhighlight> |
||
Function for performing mid-point rectangular integration: midPointRectIntegration.m |
Function for performing mid-point rectangular integration: midPointRectIntegration.m |
||
< |
<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</ |
end</syntaxhighlight> |
||
Function for performing trapezoidal integration: trapezoidalIntegration.m |
Function for performing trapezoidal integration: trapezoidalIntegration.m |
||
< |
<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</ |
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"). |
||
< |
<syntaxhighlight lang="matlab">integral = quad(f,a,b,tol)</syntaxhighlight> |
||
Using anonymous functions |
Using anonymous functions |
||
< |
<syntaxhighlight lang="matlab">trapezoidalIntegration(@(x)( exp(-(x.^2)) ),0,10,100000) |
||
ans = |
ans = |
||
0.886226925452753</ |
0.886226925452753</syntaxhighlight> |
||
Using predefined functions |
Using predefined functions |
||
Built-in MATLAB function sin(x): |
Built-in MATLAB function sin(x): |
||
< |
<syntaxhighlight lang="matlab">quad(@sin,0,pi,1/1000000000000) |
||
ans = |
ans = |
||
2.000000000000000</ |
2.000000000000000</syntaxhighlight> |
||
User defined scripts and functions: |
User defined scripts and functions: |
||
fermiDirac.m |
fermiDirac.m |
||
< |
<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</ |
end</syntaxhighlight> |
||
< |
<syntaxhighlight lang="matlab"> rightRectIntegration(@fermiDirac,-1,1,1000000) |
||
ans = |
ans = |
||
0.999998006023282</ |
0.999998006023282</syntaxhighlight> |
||
=={{header|Maxima}}== |
=={{header|Maxima}}== |
||
< |
<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;</ |
trapezium(1/x, x, 1, 100, 10000) - log(100), bfloat;</syntaxhighlight> |
||
=={{header|Nim}}== |
=={{header|Nim}}== |
||
{{trans|Python}} |
{{trans|Python}} |
||
< |
<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)</ |
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:< |
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.</ |
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. ) |
||
]</ |
]</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) |
||
]</ |
]</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</ |
) 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>. |
||
< |
<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)</ |
test(x->x, 0, 6000, 6000000)</syntaxhighlight> |
||
Results: |
Results: |
||
Line 3,880: | Line 3,880: | ||
=={{header|Pascal}}== |
=={{header|Pascal}}== |
||
< |
<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;</ |
end;</syntaxhighlight> |
||
=={{header|Perl}}== |
=={{header|Perl}}== |
||
{{trans|Raku}} |
{{trans|Raku}} |
||
< |
<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);</ |
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}}== |
||
<!--< |
<!--<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> |
||
<!--</ |
<!--</syntaxhighlight>--> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 4,127: | Line 4,127: | ||
=={{header|PicoLisp}}== |
=={{header|PicoLisp}}== |
||
< |
<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)))</ |
(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}}== |
||
< |
<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}}== |
||
< |
<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$) </ |
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. |
||
< |
<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))))</ |
float(integrate( func, a, b, steps, rule))))</syntaxhighlight> |
||
'''Tests''' |
'''Tests''' |
||
< |
<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))))</ |
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 |
||
< |
<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)</ |
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. |
||
< |
<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</ |
# 1.8e+07 1.8e+07 1.8e+07 1.8e+07 1.8e+07</syntaxhighlight> |
||
=={{header|Racket}}== |
=={{header|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: |
||
< |
<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 |
<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;</ |
.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: there was virtually no difference in accuracy between '''numeric digits 9''' (the default) and '''numeric digits 20'''. |
Note: there was virtually no difference in accuracy between '''numeric digits 9''' (the default) and '''numeric digits 20'''. |
||
< |
<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</ |
return $*h/2</syntaxhighlight> |
||
{{out|output|text= when using the default inputs:}} |
{{out|output|text= when using the default inputs:}} |
||
<pre> |
<pre> |
||
Line 4,831: | Line 4,831: | ||
=={{header|Ring}}== |
=={{header|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}} |
||
< |
<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</ |
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. |
||
< |
<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)); |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 5,031: | Line 5,031: | ||
=={{header|Scala}}== |
=={{header|Scala}}== |
||
< |
<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) |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
Output: |
Output: |
||
<pre>rectangular left : 0,245025 |
<pre>rectangular left : 0,245025 |
||
Line 5,093: | Line 5,093: | ||
=={{header|Scheme}}== |
=={{header|Scheme}}== |
||
< |
<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))</ |
(define s (integrate square 0 1 10 simpson))</syntaxhighlight> |
||
=={{header|SequenceL}}== |
=={{header|SequenceL}}== |
||
< |
<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]);</ |
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}} |
||
< |
<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);</ |
tryem('x', { _ }, 0, 6_000, 6_000_000, 18_000_000);</syntaxhighlight> |
||
=={{header|Standard ML}}== |
=={{header|Standard ML}}== |
||
< |
<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 )</ |
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</ |
end</syntaxhighlight> |
||
'''Output''' |
'''Output''' |
||
Line 5,334: | Line 5,334: | ||
=={{header|Swift}}== |
=={{header|Swift}}== |
||
< |
<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 }) }))</ |
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}}== |
||
< |
<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] |
||
} |
} |
||
}</ |
}</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>. |
||
< |
<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"</ |
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. |
||
< |
<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"</ |
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. |
||
< |
<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"></ |
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. |
||
< |
<syntaxhighlight lang="ursala">#cast %eL |
||
examples = <.left,midpoint,rignt,trapezium,simpson> (sqr,0.,pi,10)</ |
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). |
||
< |
<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</ |
End Sub</syntaxhighlight> |
||
=={{header|Wren}}== |
=={{header|Wren}}== |
||
{{trans|Kotlin}} |
{{trans|Kotlin}} |
||
{{libheader|Wren-fmt}} |
{{libheader|Wren-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}}== |
||
< |
<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); |
||
]</ |
]</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 |
||
< |
<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</ |
end sub</syntaxhighlight> |
||
=={{header|zkl}}== |
=={{header|zkl}}== |
||
{{trans|D}} |
{{trans|D}} |
||
< |
<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(); |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |