Numerical integration: Difference between revisions

This solution creates a generic package into which the function F(X) is passed during generic instantiation. The first part is the package specification. The second part is the package body.

 

with function F(X : Long_Float) return Long_Float;

package Integrate is

end Simpson;

=={{header|ALGOL 68}}==

###############

h / 6 * (f(a) + f(b) + 4 * sum1 + 2 * sum2)

END # simpson #;

 

=={{header|BASIC}}==

{{trans|Java}}

 

h = (b - a) / n

sum = 0

simpson = h / 6 * (f(a) + f(b) + 4 * sum1 + 2 * sum2)

 

=={{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)

double rr = integrate(f, 0.0, 1.0, 10, rectangular(rectangular::right));

double t = integrate(f, 0.0, 1.0, 10, trapezium());

 

=={{header|D}}==

The function f(x) is Func1/Func2 in this program.

import std.stdio ;

import std.math ;

writefln("%12.9f (%3dstep).",

Sigma(&Func2, 512, -1.0L, 2.0L, Sigma.Method.SIMP),512) ;

Parts of the output:

<pre>Function(x) = 2 / (1 + 4x^2)

 

=={{header|Forth}}==

defer method ( fn F: x -- fn[x] )

test mid fn2 \ 2.496091

test trap fn2 \ 2.351014

 

=={{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)

 

Use Array Initializers, Pointers, Array invocation of Elemental functions, Elemental array-array and array-scalar arithmetic, and the SUM intrinsic function

integer, parameter :: n = 20 ! or whatever

real, parameter :: a = 0.0, b = 15.0 ! or whatever

! Simpson's rule integral

simpson = h / 6 * sum(fleft + fright + 4*fmid)

 

=={{header|Haskell}}==

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''.

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

ws = concat $ replicate n vs

c = a + h/2

 

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

inter n [] = x : inter (n-1) xs

inter n [y] = (x+y) : inter (n-1) xs

 

And now we can just define

 

intRightRect = integrateClosed 1 [0,1]

intMidRect = integrateOpen 1 [1]

intTrapezium = integrateClosed 2 [1,1]

 

or, as easily, some additional schemes:

 

intOpen1 = integrateOpen 2 [3,3]

 

Some examples:

=={{header|Java}}==

The function in this example is assumed to be <tt>f(double x)</tt>.

public static double leftRect(double a, double b, double n){

double h = (b - a) / n;

}

//assume f(double x) elsewhere in the class

 

=={{header|Logo}}==

output invoke :fn :x

end

print integrate "i.mid "fn2 4 -1 2 ; 2.496091

print integrate "i.trapezium "fn2 4 -1 2 ; 2.351014

 

=={{header|OCaml}}==

let h = (b -. a) /. float_of_int steps in

let rec helper i s =

let rr = integrate square 0. 1. 10 right_rect

let t = integrate square 0. 1. 10 trapezium

 

=={{header|Pascal}}==

=={{header|Scheme}}==

 

(define h (/ (- b a) steps))

(* h

(define rr (integrate square 0 1 10 right-rect))

(define t (integrate square 0 1 10 trapezium))