Numerical integration/Adaptive Simpson's method: Difference between revisions
m (→{{header|C}}: Corrected minor error.) |
(Added Kotlin) |
||
Line 130: | Line 130: | ||
<pre> |
<pre> |
||
Simpson's integration of sine from 0 to 1 = 0.459698 |
Simpson's integration of sine from 0 to 1 = 0.459698 |
||
</pre> |
|||
=={{header|Kotlin}}== |
|||
{{trans|Go}} |
|||
<lang scala>// Version 1.2.71 |
|||
import kotlin.math.abs |
|||
import kotlin.math.sin |
|||
typealias F = (Double) -> Double |
|||
typealias T = Triple<Double, Double, Double> |
|||
/* "structured" adaptive version, translated from Racket */ |
|||
fun quadSimpsonsMem(f: F, a: Double, fa: Double, b: Double, fb: Double): T { |
|||
// Evaluates Simpson's Rule, also returning m and f(m) to reuse |
|||
val m = (a + b) / 2 |
|||
val fm = f(m) |
|||
val simp = abs(b - a) / 6 * (fa + 4 * fm + fb) |
|||
return T(m, fm, simp) |
|||
} |
|||
fun quadAsrRec(f: F, a: Double, fa: Double, b: Double, fb: Double, |
|||
eps: Double, whole: Double, m: Double, fm: Double): Double { |
|||
// Efficient recursive implementation of adaptive Simpson's rule. |
|||
// Function values at the start, middle, end of the intervals are retained. |
|||
val (lm, flm, left) = quadSimpsonsMem(f, a, fa, m, fm) |
|||
val (rm, frm, right) = quadSimpsonsMem(f, m, fm, b, fb) |
|||
val delta = left + right - whole |
|||
if (abs(delta) <= eps * 15) return left + right + delta / 15 |
|||
return quadAsrRec(f, a, fa, m, fm, eps / 2, left, lm, flm) + |
|||
quadAsrRec(f, m, fm, b, fb, eps / 2, right, rm, frm) |
|||
} |
|||
fun quadAsr(f: F, a: Double, b: Double, eps: Double): Double { |
|||
// Integrate f from a to b using ASR with max error of eps. |
|||
val fa = f(a) |
|||
val fb = f(b) |
|||
val (m, fm, whole) = quadSimpsonsMem(f, a, fa, b, fb) |
|||
return quadAsrRec(f, a, fa, b, fb, eps, whole, m, fm) |
|||
} |
|||
fun main(args: Array<String>) { |
|||
val a = 0.0 |
|||
val b = 1.0 |
|||
val sinx = quadAsr(::sin, a, b, 1.0e-09) |
|||
println("Simpson's integration of sine from $a to $b = ${"%6f".format(sinx)}") |
|||
}</lang> |
|||
{{output}} |
|||
<pre> |
|||
Simpson's integration of sine from 0.0 to 1.0 = 0.459698 |
|||
</pre> |
</pre> |
||
Revision as of 13:29, 26 October 2018
Lychee (1969)'s Modified Adaptive Simpson's method (doi:10.1145/321526.321537) is a numerical quadrature method that recursively bisects the interval until the precision is high enough.
; Lychee's ASR, Modifications 1, 2, 3 procedure _quad_asr_simpsons(f, a, fa, b, fb) m := (a + b) / 2 fm := f(m) h := b - a return multiple [m, fm, (h / 6) * (f(a) + f(b) + 4*sum1 + 2*sum2)] procedure _quad_asr(f, a, fa, b, fb, tol, whole, m, fm, depth) lm, flm, left := _quad_asr_simpsons(f, a, fa, m, fm) rm, frm, right := _quad_asr_simpsons(f, m, fm, b, fb) delta := left + right - whole tol' := tol / 2 if depth <= 0 or tol' == tol or abs(delta) <= 15 * tol: return left + right + delta / 15 else: return _quad_asr(f, a, fa, m, fm, tol', left , lm, flm, depth - 1) + _quad_asr(f, m, fm, b, fb, tol', right, rm, frm, depth - 1) procedure quad_asr(f, a, b, tol, depth) fa := f(a) fb := f(b) m, fm, whole := _quad_asr_simpsons(f, a, fa, b, fb) return _quad_asr(f, a, fa, b, fb, tol, whole, m, fm, depth) |
C
<lang c>#include <stdio.h>
- include <math.h>
typedef struct { double m; double fm; double simp; } triple;
/* "structured" adaptive version, translated from Racket */ triple _quad_simpsons_mem(double (*f)(double), double a, double fa, double b, double fb) {
// Evaluates Simpson's Rule, also returning m and f(m) to reuse. double m = (a + b) / 2; double fm = f(m); double simp = fabs(b - a) / 6 * (fa + 4*fm + fb); triple t = {m, fm, simp}; return t;
}
double _quad_asr(double (*f)(double), double a, double fa, double b, double fb, double eps, double whole, double m, double fm) {
// Efficient recursive implementation of adaptive Simpson's rule. // Function values at the start, middle, end of the intervals are retained. triple lt = _quad_simpsons_mem(f, a, fa, m, fm); triple rt = _quad_simpsons_mem(f, m, fm, b, fb); double delta = lt.simp + rt.simp - whole; if (fabs(delta) <= eps * 15) return lt.simp + rt.simp + delta/15; return _quad_asr(f, a, fa, m, fm, eps/2, lt.simp, lt.m, lt.fm) + _quad_asr(f, m, fm, b, fb, eps/2, rt.simp, rt.m, rt.fm);
}
double quad_asr(double (*f)(double), double a, double b, double eps) {
// Integrate f from a to b using ASR with max error of eps. double fa = f(a); double fb = f(b); triple t = _quad_simpsons_mem(f, a, fa, b, fb); return _quad_asr(f, a, fa, b, fb, eps, t.simp, t.m, t.fm);
}
int main(){
double a = 0.0, b = 1.0; double sinx = quad_asr(sin, a, b, 1e-09); printf("Simpson's integration of sine from %g to %g = %f\n", a, b, sinx); return 0;
}</lang>
- Output:
Simpson's integration of sine from 0 to 1 = 0.459698
Go
Like the zkl entry, this is also a translation of the Python code in the Wikipedia article. <lang go>package main
import (
"fmt" "math"
)
type F = func(float64) float64
/* "structured" adaptive version, translated from Racket */ func quadSimpsonsMem(f F, a, fa, b, fb float64) (m, fm, simp float64) {
// Evaluates Simpson's Rule, also returning m and f(m) to reuse. m = (a + b) / 2 fm = f(m) simp = math.Abs(b-a) / 6 * (fa + 4*fm + fb) return
}
func quadAsrRec(f F, a, fa, b, fb, eps, whole, m, fm float64) float64 {
// Efficient recursive implementation of adaptive Simpson's rule. // Function values at the start, middle, end of the intervals are retained. lm, flm, left := quadSimpsonsMem(f, a, fa, m, fm) rm, frm, right := quadSimpsonsMem(f, m, fm, b, fb) delta := left + right - whole if math.Abs(delta) <= eps*15 { return left + right + delta/15 } return quadAsrRec(f, a, fa, m, fm, eps/2, left, lm, flm) + quadAsrRec(f, m, fm, b, fb, eps/2, right, rm, frm)
}
func quadAsr(f F, a, b, eps float64) float64 {
// Integrate f from a to b using ASR with max error of eps. fa, fb := f(a), f(b) m, fm, whole := quadSimpsonsMem(f, a, fa, b, fb) return quadAsrRec(f, a, fa, b, fb, eps, whole, m, fm)
}
func main() {
a, b := 0.0, 1.0 sinx := quadAsr(math.Sin, a, b, 1e-09) fmt.Printf("Simpson's integration of sine from %g to %g = %f\n", a, b, sinx)
}</lang>
- Output:
Simpson's integration of sine from 0 to 1 = 0.459698
Kotlin
<lang scala>// Version 1.2.71
import kotlin.math.abs import kotlin.math.sin
typealias F = (Double) -> Double typealias T = Triple<Double, Double, Double>
/* "structured" adaptive version, translated from Racket */ fun quadSimpsonsMem(f: F, a: Double, fa: Double, b: Double, fb: Double): T {
// Evaluates Simpson's Rule, also returning m and f(m) to reuse val m = (a + b) / 2 val fm = f(m) val simp = abs(b - a) / 6 * (fa + 4 * fm + fb) return T(m, fm, simp)
}
fun quadAsrRec(f: F, a: Double, fa: Double, b: Double, fb: Double,
eps: Double, whole: Double, m: Double, fm: Double): Double { // Efficient recursive implementation of adaptive Simpson's rule. // Function values at the start, middle, end of the intervals are retained. val (lm, flm, left) = quadSimpsonsMem(f, a, fa, m, fm) val (rm, frm, right) = quadSimpsonsMem(f, m, fm, b, fb) val delta = left + right - whole if (abs(delta) <= eps * 15) return left + right + delta / 15 return quadAsrRec(f, a, fa, m, fm, eps / 2, left, lm, flm) + quadAsrRec(f, m, fm, b, fb, eps / 2, right, rm, frm)
}
fun quadAsr(f: F, a: Double, b: Double, eps: Double): Double {
// Integrate f from a to b using ASR with max error of eps. val fa = f(a) val fb = f(b) val (m, fm, whole) = quadSimpsonsMem(f, a, fa, b, fb) return quadAsrRec(f, a, fa, b, fb, eps, whole, m, fm)
}
fun main(args: Array<String>) {
val a = 0.0 val b = 1.0 val sinx = quadAsr(::sin, a, b, 1.0e-09) println("Simpson's integration of sine from $a to $b = ${"%6f".format(sinx)}")
}</lang>
- Output:
Simpson's integration of sine from 0.0 to 1.0 = 0.459698
zkl
<lang zkl># "structured" adaptive version, translated from Racket fcn _quad_simpsons_mem(f, a,fa, b,fb){
#Evaluates the Simpson's Rule, also returning m and f(m) to reuse""" m,fm := (a + b)/2, f(m); return(m,fm, (b - a).abs()/6*(fa + fm*4 + fb));
} fcn _quad_asr(f, a,fa, b,fb, eps, whole, m,fm){
# Efficient recursive implementation of adaptive Simpson's rule. # Function values at the start, middle, end of the intervals are retained.
lm,flm,left := _quad_simpsons_mem(f, a,fa, m,fm); rm,frm,right := _quad_simpsons_mem(f, m,fm, b,fb); delta:=left + right - whole; if(delta.abs() <= eps*15) return(left + right + delta/15); _quad_asr(f, a,fa, m,fm, eps/2, left , lm,flm) + _quad_asr(f, m,fm, b,fb, eps/2, right, rm,frm)
} fcn quad_asr(f,a,b,eps){
#Integrate f from a to b using Adaptive Simpson's Rule with max error of eps fa,fb := f(a),f(b); m,fm,whole := _quad_simpsons_mem(f, a,fa, b,fb); _quad_asr(f, a,fa, b,fb, eps,whole,m,fm);
}</lang> <lang zkl>sinx:=quad_asr((1.0).sin.unbind(), 0.0, 1.0, 1e-09); println("Simpson's integration of sine from 1 to 2 = ",sinx);</lang>
- Output:
Simpson's integration of sine from 1 to 2 = 0.459698