Skip to content

Instantly share code, notes, and snippets.

@Denommus
Last active April 22, 2022 15:24
Show Gist options
  • Save Denommus/6370332 to your computer and use it in GitHub Desktop.
Save Denommus/6370332 to your computer and use it in GitHub Desktop.
Numeric integral implementation (Simpson method) in different languages
integral :: (Fractional a, Ord a) => (a -> a) -> Integer -> a -> a -> a
integral f p a b
| a==b = 0
| otherwise = total 0 a
where dx = (b-a)/fromInteger p
total t x | x2>b = t
| otherwise = total (t+(dx*(f x+(4*f ((x+x2)/2))+f x2)/6)) x2
where x2 = x+dx
(defun integral (f precision)
(lambda (x1 x2)
(let ((dx (/ (- x2 x1) precision)))
(loop for a = x1 then (+ a dx)
for b = (+ a dx) until (> b x2)
summing (/ (* (- b a)
(+ (funcall f a)
(* 4 (funcall f (/ (+ a b) 2)))
(funcall f b)))
6) into total
finally (return total)))))
let integral f p a b =
if a==b then 0.
else let dx = (b-.a)/.float_of_int p in
let rec total t x = let x2=x+.dx in
if x2>b then t
else total (t+.(dx*.(f x+.(4.*.f ((x+.x2)/.2.))+.f x2)/.6.)) x2 in
total 0. a
def integral(f, precision)
->(x1, x2) do
dx = (x2-x1)/precision
a = x1
b = a+dx
total = 0
while b <= x2 do
total += ((b - a) *
(f.call(a) +
4 * f.call((a + b) / 2.0) +
f.call(b)))/6.0
a += dx
b += dx
end
total
end
end
// After trying lots of solutions, I ended up with something that doesn't return a closure
// For now, all the kinds of closures that Rust gives me are impractical for this use-case.
fn integral(f: &|f32|->f32, p: u32, a: f32, b: f32) -> f32 {
if p == 1 {
(b-a) * ((*f)(a) + 4.0 * (*f)((a+b)/2.0) + (*f)(b))/6.0
}
else {
let mid = (a+b)/2.0;
integral(f, p-1, a, mid) + integral(f, p-1, mid, b)
}
}
fn main() {
println(integral(&|x| x*x, 10, 1.0, 2.0).to_str());
}
(define (integral f precision)
(lambda (x1 x2)
(let ([dx (/ (- x2 x1) precision)])
(define (loop a b total)
(if (> b x2)
total
(loop (+ a dx)
(+ b dx)
(+ total
(/ (* (- b a)
(+ (f a)
(* 4 (f (/ (+ a b) 2)))
(f b)))
6)))))
(loop x1 (+ x1 dx) 0))))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment