Numerical integration — trapezoidal and Simpson's rule
Approximate the area under f over [a,b] with n panels. The trapezoid rule averages the endpoints of each panel; Simpson's rule fits parabolas and is exact for cubics. Here f(x) = x², whose integral over [0,1] is 1/3.
A Soma program from the math part of the verified corpus: 66 lines, cells Integration, IntegrationTests. It passes soma check with soma 3.2.4.
How it is built
The cell Integration declares 4 signals in its face: f(x: Float) -> Float, trapezoid(a: Float, b: Float, n: Int) -> Float, simpson(a: Float, b: Float, n: Int) -> Float and run() -> Map. Handlers: f (line 20), trapezoid (line 24), simpson (line 35) and run (line 48).
The file also holds 1 test cell (IntegrationTests) with 3 assertions, run by soma test.
The source
integration.cell
// ============================================================================
// Numerical integration — trapezoidal and Simpson's rule
// ============================================================================
//
// Approximate the area under f over [a,b] with n panels. The trapezoid rule
// averages the endpoints of each panel; Simpson's rule fits parabolas and is
// exact for cubics. Here f(x) = x², whose integral over [0,1] is 1/3.
//
// soma run integration.cell run
// ============================================================================
cell Integration {
face {
signal f(x: Float) -> Float
signal trapezoid(a: Float, b: Float, n: Int) -> Float
signal simpson(a: Float, b: Float, n: Int) -> Float
signal run() -> Map
}
on f(x: Float) [native] {
return x * x
}
on trapezoid(a: Float, b: Float, n: Int) [native] {
let h = (b - a) / to_float(n)
let s = (f(a) + f(b)) / 2.0
let i = 1
while i < n {
s = s + f(a + to_float(i) * h)
i = i + 1
}
return s * h
}
on simpson(a: Float, b: Float, n: Int) [native] {
// n must be even
let h = (b - a) / to_float(n)
let s = f(a) + f(b)
let i = 1
while i < n {
let x = a + to_float(i) * h
if i % 2 == 1 { s = s + 4.0 * f(x) } else { s = s + 2.0 * f(x) }
i = i + 1
}
return s * h / 3.0
}
on run() {
return map(
"trapezoid", trapezoid(0.0, 1.0, 1000),
"simpson", simpson(0.0, 1.0, 100)
)
}
}
cell test IntegrationTests {
rules {
// integral of x^2 over [0,1] is 1/3
assert abs(trapezoid(0.0, 1.0, 1000) - 0.3333333333333) < 0.0001
// Simpson is exact for polynomials up to cubic
assert abs(simpson(0.0, 1.0, 100) - 0.3333333333333) < 0.0000001
// integral over [0,3] of x^2 = 9
assert abs(simpson(0.0, 3.0, 100) - 9.0) < 0.0000001
}
}What the toolchain says
$ soma check integration.cell
✓ All checks passed.$ soma test integration.cell
[native] compiling 3 handler(s) for cell 'Integration'...
[native] compiled → .soma_cache/native/native_integration_499e01a06da825f4.dylib
test IntegrationTests ...
✓ assert abs(trapezoid(0.0, 1.0, 1000) - 0.3333333333333) < 0.0001
✓ assert abs(simpson(0.0, 1.0, 100) - 0.3333333333333) < 0.0000001
✓ assert abs(simpson(0.0, 3.0, 100) - 9.0) < 0.0000001
3 tests: 3 passed, 0 failed$ soma verify integration.cell
Verifying integration.cell...
Integration/termination:
⚠ handler `trapezoid`: while-loop without provable termination bound (write `while [loop_bound(N)] cond { … }`, or a for-loop over a bounded range)
⚠ handler `simpson`: while-loop without provable termination bound (write `while [loop_bound(N)] cond { … }`, or a for-loop over a bounded range)
⚠ `run` calls `Integration.simpson`, which may not terminate
0 passed, 3 warnings, 0 failures
Temporal: 0 passed, 0 failed
VERIFY OK$ soma run integration.cell run
[native] using cached native_integration_499e01a06da825f4.dylib for cell 'Integration'
{"trapezoid": 0.33333349999999995, "simpson": 0.3333333333333334}