Numerical integration: Difference between revisions
m
syntax highlighting fixup automation
m (→{{header|R}}) |
Thundergnat (talk | contribs) m (syntax highlighting fixup automation) |
||
Line 48:
{{trans|Nim}}
<
R f(x)
Line 87:
(simpson, ‘simpson’)]
print("#. integrated using #.\n from #. to #. (#. steps) = #.".format(
func_name, rule_name, a, b, steps, integrate(func, a, b, steps, rule)))</
{{out}}
Line 135:
=={{header|ActionScript}}==
Integration functions:
<
{
var sum:Number = 0;
Line 185:
}
return (dx/6) * (f(a) + f(b) + 4*sum1 + 2*sum2);
}</
Usage:
<
return (2/(1+ 4*(n*n)));
}
Line 195:
trace(trapezium(f1, -1, 2 ,4 ));
trace(simpson(f1, -1, 2 ,4 ));
</syntaxhighlight>
=={{header|Ada}}==
Specification of a generic package implementing the five specified kinds of numerical integration:
<
type Scalar is digits <>;
with function F (X : Scalar) return Scalar;
Line 208:
function Trapezium (A, B : Scalar; N : Positive) return Scalar;
function Simpsons (A, B : Scalar; N : Positive) return Scalar;
end Integrate;</
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:
<
function Left_Rectangular (A, B : Scalar; N : Positive) return Scalar is
H : constant Scalar := (B - A) / Scalar (N);
Line 275:
return (H / 3.0) * (F (A) + F (B) + 4.0 * Sum_U + 2.0 * Sum_E);
end Simpsons;
end Integrate;</
Test driver:
<
with Integrate;
Line 382:
end X;
end Numerical_Integration;
</syntaxhighlight>
=={{header|ALGOL 68}}==
<
###############
Line 502:
test integrators( "x ", ( LONG REAL x )LONG REAL: x, 0, 6 000, 6 000 000 );
SKIP</
{{out}}
<pre>
Line 513:
=={{header|ALGOL W}}==
{{Trans|ALGOL 68}}
<
long real procedure leftRect ( long real procedure f
Line 641:
testIntegrators2( "x ", xValue, 0, 6000, 6000000 )
end
end.</
=={{header|AutoHotkey}}==
ahk [http://www.autohotkey.com/forum/viewtopic.php?t=44657&postdays=0&postorder=asc&start=139 discussion]
<
MsgBox % Rect("fun", 0, 1, 10) ; 0.50 mid
MsgBox % Rect("fun", 0, 1, 10, 1) ; 0.55 right
Line 679:
fun(x) { ; linear test function
Return x
}</
=={{header|BASIC}}==
{{works with|QuickBasic|4.5}}
{{trans|Java}}
<
h = (b - a) / n
sum = 0
Line 734:
simpson = h / 6 * (f(a) + f(b) + 4 * sum1 + 2 * sum2)
END FUNCTION</
=={{header|BBC BASIC}}==
<
@% = 12 : REM Column width
Line 807:
NEXT
x = a
= (d / 6) * (f + EVAL(x$) + 4 * s1 + 2 * s2)</
'''Output:'''
<pre>
Line 818:
=={{header|C}}==
<
#include <stdlib.h>
#include <math.h>
Line 875:
return h / 6.0 * (func(from) + func(to) + 4.0 * sum1 + 2.0 * sum2);
}</
<
double f3(double x)
{
Line 947:
printf("\n");
}
}</
=={{header|C sharp|C#}}==
<
using System.Collections.Generic;
using System.Linq;
Line 1,068:
TestApproximationMethods(new DefiniteIntegral(x => x, new Interval(0, 6000)), 6000000);
}
}</
Output:
<syntaxhighlight lang="text">0.2499500025
0.24999999875
0.2500500025
Line 1,089:
18000003
18000000
18000000</
=={{header|C++}}==
Due to their similarity, it makes sense to make the integration method a policy.
<
template<typename Method, typename F, typename Float>
double integrate(F f, Float a, Float b, int steps, Method m)
Line 1,156:
double rr = integrate(f, 0.0, 1.0, 10, rectangular(rectangular::right));
double t = integrate(f, 0.0, 1.0, 10, trapezium());
double s = integrate(f, 0.0, 1.0, 10, simpson());</
=={{header|Chapel}}==
<
proc f1(x:real):real {
return x**3;
Line 1,289:
writeln("simpsonsIntegration: calculated = ", calculated, "; exact = ", exact, "; difference = ", abs(calculated - exact));
writeln();
</syntaxhighlight>
output
<syntaxhighlight lang="text">
f(x) = x**3 with 100 steps from 0 to 1
leftRectangleIntegration: calculated = 0.245025; exact = 0.25; difference = 0.004975
Line 1,319:
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
</syntaxhighlight>
=={{header|CoffeeScript}}==
{{trans|python}}
<
rules =
left_rect: (f, x, h) -> f(x)
Line 1,357:
result = integrate func, a, b, steps, rule
console.log rule_name, result
</syntaxhighlight>
output
<syntaxhighlight lang="text">
> coffee numerical_integration.coffee
-- tests for cube with 100 steps from 0 to 1
Line 1,385:
trapezium 17999999.999999993
simpson 17999999.999999993
</syntaxhighlight>
=={{header|Comal}}==
{{works with|OpenComal on Linux}}
<syntaxhighlight lang="comal">
1000 PRINT "F(X)";" FROM";" TO";" L-Rect";" M-Rect";" R-Rect ";" Trapez";" Simpson"
1010 fromval:=0
Line 1,484:
1920 RETURN x
1930 ENDFUNC
</syntaxhighlight>
{{out}}
<pre>
Line 1,496:
=={{header|Common Lisp}}==
<
(* d (loop for x from a below b by d summing (funcall f x))))
Line 1,522:
(funcall f b)
(* 4 sum1)
(* 2 sum2))))))</
=={{header|D}}==
<
template integrate(alias method) {
Line 1,587:
writeln();
}
}</
Output:
<pre>rectangular left: 0.202500
Line 1,614:
===A faster version===
This version avoids function pointers and delegates, same output:
<
template integrate(alias method) {
Line 1,684:
writeln();
}
}</
=={{header|Delphi}}==
{{libheader| System.SysUtils}}
{{Trans|Python}}
<
{$APPTYPE CONSOLE}
Line 1,785:
end;
readln;
end.</
{{out}}
<pre>Integrate x^3
Line 1,815:
{{trans|Python}}
<
def leftRect(f, x, h) {
Line 1,840:
def h := (b-a) / steps
return h * accum 0 for i in 0..!steps { _ + meth(f, a+i*h, h) }
}</
<
# value: 105.33333333333334
? integrate(fn x { x ** 9 }, 0, 1, 300, simpson)
# value: 0.10000000002160479</
=={{header|Elixir}}==
<
@funs ~w(leftrect midrect rightrect trapezium simpson)a
Line 1,884:
f4 = fn x -> x end
IO.puts "\nf(x) = x, where x is [0,6000], with 6,000,000 approximations."
Numerical.integrate(f4, 0, 6000, 6_000_000)</
{{out}}
Line 1,918:
=={{header|Euphoria}}==
<
atom h, sum
h = (bounds[2]-bounds[1])/n
Line 1,996:
? int_rightrect({0,10},1000,routine_id("x"))
? int_midrect({0,10},1000,routine_id("x"))
? int_simpson({0,10},1000,routine_id("x"))</
Output:
Line 2,016:
=={{header|F Sharp}}==
<
// integration methods
let left f dx x = f x * dx
Line 2,042:
|> Seq.map (fun ((f, a, b, n), method) -> integrate a b f n method)
|> Seq.iter (printfn "%f")
</syntaxhighlight>
=={{header|Factor}}==
<
USE: math.functions
IN: scratchpad 0 1 [ 3 ^ ] integrate-simpson .
Line 2,057:
IN: scratchpad 6000000 num-steps set-global
IN: scratchpad 0 6000 [ ] integrate-simpson .
18000000</
=={{header|Forth}}==
<
defer method ( fn F: x -- fn[x] )
Line 2,099:
test mid fn2 \ 2.496091
test trap fn2 \ 2.351014
test simpson fn2 \ 2.447732</
=={{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:
<
real :: elemf, x
elemf = f(x)
end function elemf</
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.
<
implicit none
Line 2,183:
end function integrate
end module Integration</
Usage example:
<
use Integration
use FunctionHolder
Line 2,197:
print *, integrate(afun, 0., 3**(1/3.), method='trapezoid')
end program IntegrationTest</
The FunctionHolder module:
<
implicit none
Line 2,213:
end function afun
end module FunctionHolder</
=={{header|FreeBASIC}}==
Based on the BASIC entry and the BBC BASIC entry
<
' compile with: fbc -s console
Line 2,341:
Print : Print "hit any key to end program"
Sleep
End</
{{out}}
<pre>function range steps leftrect midrect rightrect trap simpson
Line 2,350:
=={{header|Go}}==
<
import (
Line 2,513:
fmt.Println("")
}
}</
{{out}}
<pre>
Line 2,554:
=={{header|Groovy}}==
Solution:
<
assert (bounds.size() == 2) && (bounds[0] instanceof Double) && (bounds[1] instanceof Double) && (nRect > 0)
}
Line 2,597:
h/3*((fLeft + fRight).sum() + 4*(fMid.sum()))
}
}</
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).
<
double sinIntegralCalculated = -(Math.cos(Math.PI) - Math.cos(0d))
Line 2,641:
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))
assert ((simpsonsIntegral([-Math.E, Math.PI], 1, cubicPoly) - cpIntegralCalcMinusEToPI)/ cpIntegralCalcMinusEToPI).abs() < tolerance**2.5 // 1 in 10 billion</
Requested Demonstrations:
<
println ([" LeftRect": leftRectIntegral([0d, 1d], 100) { it**3 }])
println (["RightRect": rightRectIntegral([0d, 1d], 100) { it**3 }])
Line 2,674:
println (["Trapezoid": trapezoidIntegral([0d, 6000d], 6000000) { it }])
println ([" Simpsons": simpsonsIntegral([0d, 6000d], 6000000) { it }])
println ()</
Output:
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
<
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:
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:
<
integrateOpen v vs f a b n = approx f xs ws * h / v where
m = fromIntegral (length vs) * n
Line 2,721:
ws = concat $ replicate n vs
c = a + h/2
xs = [c + h * fromIntegral i | i <- [0..m-1]]</
Similarly for the closed formulas, but we need an additional function ''overlap'' which sums the coefficients overlapping at the interior interval boundaries:
<
integrateClosed v vs f a b n = approx f xs ws * h / v where
m = fromIntegral (length vs - 1) * n
Line 2,738:
inter n [] = x : inter (n-1) xs
inter n [y] = (x+y) : inter (n-1) xs
inter n (y:ys) = y : inter n ys</
And now we can just define
<
intRightRect = integrateClosed 1 [0,1]
intMidRect = integrateOpen 1 [1]
intTrapezium = integrateClosed 2 [1,1]
intSimpson = integrateClosed 3 [1,4,1]</
or, as easily, some additional schemes:
<
intOpen1 = integrateOpen 2 [3,3]
intOpen2 = integrateOpen 3 [8,-4,8]</
Some examples:
Line 2,769:
The whole program:
<
:: Fractional a
=> (a1 -> a) -> [a1] -> [a] -> a
Line 2,857:
integrations
where
indent n = take n . (++ replicate n ' ')</
{{Out}}
<pre>f(x) = x^3 [0.0,1.0] (100 approximations)
Line 2,890:
=={{header|J}}==
===Solution:===
<
'a b steps'=. 3{.y,128
size=. (b - a)%steps
Line 2,900:
trapezium=: adverb def '-: +/ u y'
simpson =: adverb def '6 %~ +/ 1 1 4 * u y, -:+/y'</
===Example usage===
====Required Examples====
<
It=: trapezium integrate
Is=: simpson integrate
Line 2,930:
1.8e7
] Is 0 6000 6e6
1.8e7</
====Older Examples====
Integrate <code>square</code> (<code>*:</code>) from 0 to π in 10 steps using various methods.
<
10.3095869962
*: trapezium integrate 0 1p1 10
10.3871026879
*: simpson integrate 0 1p1 10
10.3354255601</
Integrate <code>sin</code> from 0 to π in 10 steps using various methods.
<
sin rectangle integrate 0 1p1 10
2.00824840791
Line 2,946:
1.98352353751
sin simpson integrate 0 1p1 10
2.00000678444</
===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:
<
0 0 0 0.333333333333
0 p.. 0 0 1x NB. or using rationals
0 0 0 1r3</
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 π :
<
10.3354255601</
That said, J also has <code>d.</code> which can [http://www.jsoftware.com/help/dictionary/dddot.htm integrate] suitable functions.
<
10.3354</
=={{header|Java}}==
<
{
Line 3,087:
}
}
</syntaxhighlight>
=={{header|Julia}}==
{{works with|Julia|0.6}}
<
h = (b - a) / n
s = f(a + h / 2)
Line 3,107:
simpson(x -> x, 0, 6000, 6_000_000)
@show rst</
{{out}}
Line 3,113:
=={{header|Kotlin}}==
<
typealias Func = (Double) -> Double
Line 3,138:
integrate(0.0, 5000.0, 5_000_000) { it }
integrate(0.0, 6000.0, 6_000_000) { it }
}</
{{out}}
Line 3,170:
Following Python's presentation
<syntaxhighlight lang="scheme">
1) FUNCTIONS
Line 3,258:
trapezium 18006000
simpson 18006000
</syntaxhighlight>
=={{header|Liberty BASIC}}==
Running the big loop value would take a VERY long time & seems unnecessary.<
while 1
read x$
Line 3,379:
end
</syntaxhighlight>
Numerical integration
Line 3,415:
=={{header|Logo}}==
<
output invoke :fn :x
end
Line 3,450:
print integrate "i.mid "fn2 4 -1 2 ; 2.496091
print integrate "i.trapezium "fn2 4 -1 2 ; 2.351014
print integrate "i.simpsons "fn2 4 -1 2 ; 2.447732</
=={{header|Lua}}==
<
local h = (b - a) / n
local x = a
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, 6000, 6000000 ) )
end</
=={{header|Mathematica}}/{{header|Wolfram Language}}==
<
Module[{sum = 0, dx = (b - a)/N, x = a, n = N} ,
For[n = N, n > 0, n--, x += dx; sum += f[x];];
Line 3,553:
For[n = 1, n < N, n++, sum1 += f[a + dx*n + dx/2];
sum2 += f[a + dx*n];];
Return [(dx/6)*(f[a] + f[b] + 4*sum1 + 2*sum2)]]</
<pre>f[x_] := x^3
g[x_] := 1/x
Line 3,574:
Function for performing left rectangular integration: leftRectIntegration.m
<
format long;
Line 3,581:
integral = width * sum( f(x(1:n-1)) );
end</
Function for performing right rectangular integration: rightRectIntegration.m
<
format long;
Line 3,591:
integral = width * sum( f(x(2:n)) );
end</
Function for performing mid-point rectangular integration: midPointRectIntegration.m
<
format long;
Line 3,601:
integral = width * sum( f( (x(1:n-1)+x(2:n))/2 ) );
end</
Function for performing trapezoidal integration: trapezoidalIntegration.m
<
format long;
Line 3,610:
integral = trapz( x,f(x) );
end</
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").
<
Using anonymous functions
<
ans =
0.886226925452753</
Using predefined functions
Built-in MATLAB function sin(x):
<
ans =
2.000000000000000</
User defined scripts and functions:
fermiDirac.m
<
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
end</
<
ans =
0.999998006023282</
=={{header|Maxima}}==
<
for i from 1 thru n do s: s + subst(x = a + i * h, e),
s * h)$
Line 3,672:
2 * log(2) - 1 - %, bfloat;
trapezium(1/x, x, 1, 100, 10000) - log(100), bfloat;</
=={{header|Nim}}==
{{trans|Python}}
<
type Rule = proc(f: Function; x, h: float): float
Line 3,720:
echo fName, " integrated using ", rName
echo " from ", a, " to ", b, " (", steps, " steps) = ",
integrate(fun, float(a), float(b), steps, rule)</
{{out}}
Line 3,766:
=={{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.
First define the integration function:<
let h = (b -. a) /. float_of_int steps in
let rec helper i s =
Line 3,772:
else helper (succ i) (s +. meth f (a +. h *. float_of_int i) h)
in
h *. helper 0 0.</
( "rect_l", fun f x _ -> f x);
( "rect_m", fun f x h -> f (x +. h /. 2.) );
Line 3,778:
( "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. )
]</
( "cubic", (fun x -> x*.x*.x), 0.0, 1.0, 100);
( "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 6e3", (fun x -> x), 0.0, 6000.0, 6_000_000)
]</
List.iter (fun (s,f,lo,hi,n) ->
Printf.printf "Testing function %s:\n" s;
Line 3,789:
Printf.printf " method %s gives %.15g\n" name (integrate f lo hi n meth)
) methods
) functions</
<pre>
Testing function cubic:
Line 3,819:
=={{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>.
<
sum(i=0,n-1,f(a+(b-a)*i/n), 0.)*(b-a)/n
};
Line 3,847:
test(x->1/x, 1, 100, 1000)
test(x->x, 0, 5000, 5000000)
test(x->x, 0, 6000, 6000000)</
Results:
Line 3,880:
=={{header|Pascal}}==
<
begin
RectLeft := f(xl)
Line 3,920:
end;
integrate := integral
end;</
=={{header|Perl}}==
{{trans|Raku}}
<
sub leftrect {
Line 4,004:
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, 6_000, 6_000_000, 18_000_000);</
{{out}}
<pre>$_ ** 3
Line 4,043:
=={{header|Phix}}==
<!--<
<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>
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>
<!--</
{{out}}
<pre>
Line 4,127:
=={{header|PicoLisp}}==
<
(de leftRect (Fun X)
Line 4,158:
(*/ H Sum 1.0) ) )
(prinl (round (integrate square 3.0 7.0 30 simpson)))</
Output:
<pre>105.333</pre>
=={{header|PL/I}}==
<
f: procedure (x, function) returns (float(18));
Line 4,234:
end integrals;
</syntaxhighlight>
<pre>
Rectangle-left Rectangle-mid Rectangle-right Trapezoid Simpson
Line 4,245:
=={{header|PureBasic}}==
<
Procedure.d LeftIntegral(Start, Stop, Steps, *func.TestFunction)
Line 4,346:
Answer$+"Trapezium="+StrD(Trapezium (0,6000,6000000,@Test3()))+#CRLF$
Answer$+"Simpson ="+StrD(Simpson (0,6000,6000000,@Test3()))
MessageRequester("Answer should be 18,000,000",Answer$) </
<pre>Left =0.2353220100
Mid =0.2401367513
Line 4,373:
=={{header|Python}}==
Answers are first given using floating point arithmatic, then using fractions, only converted to floating point on output.
<
def left_rect(f,x,h):
Line 4,427:
print('%s integrated using %s\n from %r to %r (%i steps and fractions) = %r' %
(func.__name__, rule.__name__, a, b, steps,
float(integrate( func, a, b, steps, rule))))</
'''Tests'''
<
for rule in (left_rect, mid_rect, right_rect, trapezium, simpson):
print('%s integrated using %s\n from %r to %r (%i steps) = %r' %
Line 4,452:
print('%s integrated using %s\n from %r to %r (%i steps and fractions) = %r' %
(func.__name__, rule.__name__, a, b, steps,
float(integrate( func, a, b, steps, rule))))</
'''Sample test Output'''
Line 4,537:
A faster Simpson's rule integrator is
<
h = (b-a)/float(steps)
a1 = a+h/2
s1 = sum( f(a1+i*h) for i in range(0,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)</
=={{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.
<
h <- (b - a) / n
s <- 0
Line 4,578:
test(\(x) x, 0, 6000, 6e6)
# rect.left rect.right rect.mid trapezoidal simpson
# 1.8e+07 1.8e+07 1.8e+07 1.8e+07 1.8e+07</
=={{header|Racket}}==
<
#lang racket
(define (integrate f a b steps meth)
Line 4,605:
(test (λ(x) x) 0. 5000. 5000000 "IDENTITY")
(test (λ(x) x) 0. 6000. 6000000 "IDENTITY")
</syntaxhighlight>
Output:
<
CUBED
left-rect: 0.24502500000000005
Line 4,635:
trapezium: 17999999.999999993
simpson: 17999999.999999993
</syntaxhighlight>
=={{header|Raku}}==
Line 4,644:
{{works with|Rakudo|2018.09}}
<syntaxhighlight lang="raku"
sub leftrect(&f, $a, $b, $n) {
Line 4,711:
.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, 6_000, 6_000_000, 18_000_000;</
{{out}}
<pre>{ $_ ** 3 }
Line 4,751:
=={{header|REXX}}==
Note: there was virtually no difference in accuracy between '''numeric digits 9''' (the default) and '''numeric digits 20'''.
<
numeric digits 20 /*use twenty decimal digits precision. */
Line 4,798:
do x=a by h for #; $= $ + (f(x) + f(x+h))
end /*x*/
return $*h/2</
{{out|output|text= when using the default inputs:}}
<pre>
Line 4,831:
=={{header|Ring}}==
<
# Project : Numerical integration
Line 4,917:
eval("result = " + x2)
return (d / 6) * (f + result + 4 * s1 + 2 * s)
</syntaxhighlight>
Output:
<pre>
Line 4,929:
=={{header|Ruby}}==
{{trans|Tcl}}
<
f.call(left)
end
Line 4,986:
printf " %-10s %s\t(%.1f%%)\n", method, int, diff
end
end</
outputs
<pre>integral of #<Method: Object#square> from 0 to 3.14159265358979 in 10 steps
Line 5,003:
=={{header|Rust}}==
This is a partial solution and only implements trapezium integration.
<
where F: Fn(f64) -> f64
{
Line 5,022:
println!("{}", integral(|x| x, 0.0..5000.0, 5_000_000));
println!("{}", integral(|x| x, 0.0..6000.0, 6_000_000));
}</
{{out}}
Line 5,031:
=={{header|Scala}}==
<
def leftRect(f:Double=>Double, a:Double, b:Double)=f(a)
def midRect(f:Double=>Double, a:Double, b:Double)=f((a+b)/2)
Line 5,065:
print(fn3, 0, 6000, 6000000)
}
}</
Output:
<pre>rectangular left : 0,245025
Line 5,093:
=={{header|Scheme}}==
<
(define h (/ (- b a) steps))
(* h
Line 5,113:
(define rr (integrate square 0 1 10 right-rect))
(define t (integrate square 0 1 10 trapezium))
(define s (integrate square 0 1 10 simpson))</
=={{header|SequenceL}}==
<
import <Utilities/Sequence.sl>;
Line 5,175:
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]);</
{{out}}
Line 5,188:
=={{header|Sidef}}==
{{trans|Raku}}
<
var s = 0;
RangeNum(start, to, from-start).each { |i|
Line 5,241:
tryem('1/x', { 1 / _ }, 1, 100, 1000, log(100));
tryem('x', { _ }, 0, 5_000, 5_000_000, 12_500_000);
tryem('x', { _ }, 0, 6_000, 6_000_000, 18_000_000);</
=={{header|Standard ML}}==
<
val h = (b - a) / real steps
fun helper (i, s) =
Line 5,266:
val rr = integrate (square, 0.0, 1.0, 10, right_rect)
val t = integrate (square, 0.0, 1.0, 10, trapezium )
val s = integrate (square, 0.0, 1.0, 10, simpson )</
=={{header|Stata}}==
<syntaxhighlight lang="text">mata
function integrate(f,a,b,n,u,v) {
s = 0
Line 5,309:
test(&id(),0,5000,5000000)
test(&id(),0,6000,6000000)
end</
'''Output'''
Line 5,334:
=={{header|Swift}}==
<
case rectangularLeft
case rectangularRight
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) = 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 }) }))</
{{out}}
Line 5,447:
=={{header|Tcl}}==
<
proc leftrect {f left right} {
Line 5,500:
puts [format " %-10s %s\t(%.1f%%)" $method $int $diff]
}
}</
<pre>integral of square(x) from 0 to 3.141592653589793 in 10 steps
leftrect 8.836788853885448 (-14.5%)
Line 5,530:
the integrand <math>f</math>, the bounds <math>(a,b)</math>, and the number of intervals <math>n</math>.
<
#import nat
#import flo
Line 5,536:
(integral_by "m") ("f","a","b","n") =
iprod ^(* ! div\float"n" minus/"b" "a",~&) ("m" "f")*ytp (ari successor "n")/"a" "b"</
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.
<
iprod ^(* ! div\float"n" minus/"b" "a",~&) ^H(*+ "m"+ -:"f"+ * ^/~& "f",~&ytp) (ari successor "n")/"a" "b"</
As mentioned in the Haskell solution, the latter choice is preferable if evaluating the integrand
is expensive.
An integrating function is defined for each method as follows.
<
right = integral_by "f". ("l","r"). "f" "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"
simpson = integral_by "f". ("l","r"). div\6. plus:-0. <"f" "l",times/4. "f" div\2. plus/"l" "r","f" "r"></
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
Line 5,555:
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.
<
examples = <.left,midpoint,rignt,trapezium,simpson> (sqr,0.,pi,10)</
output:
<pre>
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).
<
Option Base 1
Line 5,628:
Next j
Next i
End Sub</
=={{header|Wren}}==
{{trans|Kotlin}}
{{libheader|Wren-fmt}}
<
var integrate = Fn.new { |a, b, n, f|
Line 5,655:
integrate.call(0, 5000, 5000000) { |v| v }
integrate.call(0, 6000, 6000000) { |v| v }
</syntaxhighlight>
{{out}}
Line 5,685:
=={{header|XPL0}}==
<
func real Func(FN, X); \Return F(X) for function number FN
Line 5,731:
Integrate(0.0, 5000.0, 3, 5_000_000);
Integrate(0.0, 6000.0, 3, 6_000_000);
]</
Interestingly, the small rounding errors creep in when millions of
Line 5,747:
=={{header|Yabasic}}==
Based on the XPL0entry and the Free BASIC entry
<
// by Jjuanhdez, 06/2022
Line 5,819:
print " ";
print str$(Area, frmt$)
end sub</
=={{header|zkl}}==
{{trans|D}}
<
h:=(b - a) / steps;
h*(0).reduce(steps,'wrap(s,i){ F(f, h*i + a, h) + s },0.0);
Line 5,846:
"%s %f".fmt(nm,integrate(f,a.xplode())).println() }, fs);
println();
}</
{{out}}
<pre>
|