Add cumtrapz/cumtrapz!

Project.toml

@@ -1,7 +1,7 @@
 name = "FLOWMath"
 uuid = "6cb5d3fb-0fe8-4cc2-bd89-9fe0b19a99d3"
 authors = ["Andrew Ning <[email protected]>"]
-version = "0.4.0"
+version = "0.4.1"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

docs/src/index.md

@@ -20,6 +20,31 @@ z = trapz(x, y)
 trapz
 ```
 
+There is also `cumtrapz` which returns the cumulative integral of `y` with respect to `x`
+
+```@example cumtrapz
+using FLOWMath: cumtrapz
+x = range(0.0, stop=pi+1e-15, step=pi/100)
+y = sin.(x)
+z = cumtrapz(x, y)
+extrema(z .- ((-cos.(x)) .- (-cos(x[1]))))  # compare to the exact answer -cos(x) - -cos(0)
+```
+
+and a version `cumtrapz!` that writes the result to the first argument
+
+```@example cumtrapz!
+using FLOWMath: cumtrapz!
+x = range(0.0, stop=pi+1e-15, step=pi/100)
+y = sin.(x)
+z = similar(y)
+cumtrapz!(z, x, y)
+extrema(z .- ((-cos.(x)) .- (-cos(x[1]))))  # compare to the exact answer -cos(x) - -cos(0)
+```
+
+```@docs
+cumtrapz
+cumtrapz!
+```
 ## Root Finding
 
 ### Brent's Method (1D functions)

src/FLOWMath.jl

@@ -7,7 +7,7 @@ export dot_cs_safe
 export atan_cs_safe
 
 include("quadrature.jl")
-export trapz
+export trapz, cumtrapz, cumtrapz!
 
 include("roots.jl")
 export brent

src/quadrature.jl

@@ -12,3 +12,26 @@ function trapz(x, y)
     end
     return integral
 end
+
+
+"""
+    cumtrapz(x, y)
+
+Cumulatively integrate `y` w.r.t `x` using the trapezoidal method, returning an array the same size as `x` and `y`.
+"""
+cumtrapz(x, y) = cumtrapz!(similar(y), x, y)
+
+
+"""
+    cumtrapz!(integral, x, y)
+
+Cumulatively integrate `y` w.r.t `x` using the trapezoidal method, writing the result to `integral`.
+"""
+function cumtrapz!(integral, x, y)
+    integral[begin] = 0
+    for i in eachindex(x, y, integral)[begin+1:end]
+        integral[i] = integral[i-1] + (x[i] - x[i-1])*0.5*(y[i] + y[i-1])
+    end
+
+    return integral
+end

test/runtests.jl

@@ -101,7 +101,41 @@ y = sin.(x)
 z = trapz(x, y)
 @test isapprox(z, 1.999835503887444, atol=1e-15)
 
-# -------------------------
+# ------ cumtrapz --------
+# tests from the matlab cumtrapz docs
+# https://www.mathworks.com/help/matlab/ref/cumtrapz.html
+y = Float64[1, 4, 9, 16, 25]
+x = 1:length(y)
+@test all(cumtrapz(x, y) .≈ [0.0, 2.5, 9.0, 21.5, 42.0])
+out = similar(y)
+cumtrapz!(out, x, y)
+@test all(cumtrapz(x, y) .≈ out)
+
+x = range(0, pi; length=6)
+y = sin.(x)
+@test all(isapprox.(cumtrapz(x, y), [0, 0.1847, 0.6681, 1.2657, 1.7491, 1.9338]; atol=1e-4))
+out = similar(y)
+cumtrapz!(out, x, y)
+@test all(cumtrapz(x, y) .≈ out)
+
+x = [1, 2.5, 7, 10]
+y = [5.2, 7.7,  9.6, 13.2]
+@test all(cumtrapz(x, y) .≈ [0, 9.6750, 48.6000, 82.8000])
+out = similar(y)
+cumtrapz!(out, x, y)
+@test all(cumtrapz(x, y) .≈ out)
+
+y = [4.8, 7.0, 10.5, 14.5]
+@test all(cumtrapz(x, y) .≈ [0, 8.8500, 48.2250, 85.7250])
+out = similar(y)
+cumtrapz!(out, x, y)
+@test all(cumtrapz(x, y) .≈ out)
+
+y = [4.9, 6.5, 10.2, 13.8]
+@test all(cumtrapz(x, y) .≈ [0, 8.5500, 46.1250, 82.1250])
+out = similar(y)
+cumtrapz!(out, x, y)
+@test all(cumtrapz(x, y) .≈ out)
 
 # ------ Brent's method ------