Add first try at cumtrapz/cumtrapz!

docs/Manifest.toml

@@ -1,5 +1,12 @@
 # This file is machine-generated - editing it directly is not advised
 
+[[ArgTools]]
+uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
+version = "1.1.1"
+
+[[Artifacts]]
+uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
+
 [[Base64]]
 uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 
@@ -21,6 +28,11 @@ git-tree-sha1 = "ed2c4abadf84c53d9e58510b5fc48912c2336fbb"
 uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
 version = "2.2.0"
 
+[[CompilerSupportLibraries_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
+version = "1.0.5+0"
+
 [[Conda]]
 deps = ["JSON", "VersionParsing"]
 git-tree-sha1 = "9a11d428dcdc425072af4aea19ab1e8c3e01c032"
@@ -39,7 +51,9 @@ uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
 
 [[DelimitedFiles]]
 deps = ["Mmap"]
+git-tree-sha1 = "9e2f36d3c96a820c678f2f1f1782582fcf685bae"
 uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+version = "1.9.1"
 
 [[Distributed]]
 deps = ["Random", "Serialization", "Sockets"]
@@ -57,6 +71,20 @@ git-tree-sha1 = "d497bcc45bb98a1fbe19445a774cfafeabc6c6df"
 uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 version = "0.24.5"
 
+[[Downloads]]
+deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"]
+uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
+version = "1.6.0"
+
+[[FLOWMath]]
+deps = ["LinearAlgebra", "OffsetArrays"]
+path = ".."
+uuid = "6cb5d3fb-0fe8-4cc2-bd89-9fe0b19a99d3"
+version = "0.3.3"
+
+[[FileWatching]]
+uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"
+
 [[FixedPointNumbers]]
 git-tree-sha1 = "4aaea64dd0c30ad79037084f8ca2b94348e65eaa"
 uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
@@ -78,14 +106,30 @@ git-tree-sha1 = "7ab9b8788cfab2bdde22adf9004bda7ad9954b6c"
 uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 version = "1.0.3"
 
+[[LibCURL]]
+deps = ["LibCURL_jll", "MozillaCACerts_jll"]
+uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"
+version = "0.6.4"
+
+[[LibCURL_jll]]
+deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
+uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"
+version = "8.4.0+0"
+
 [[LibGit2]]
+deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
 uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
 
+[[LibSSH2_jll]]
+deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
+uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"
+version = "1.11.0+1"
+
 [[Libdl]]
 uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
 
 [[LinearAlgebra]]
-deps = ["Libdl"]
+deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"]
 uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 
 [[Logging]]
@@ -101,9 +145,38 @@ version = "0.5.3"
 deps = ["Base64"]
 uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
 
+[[MbedTLS_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
+version = "2.28.2+0"
+
 [[Mmap]]
 uuid = "a63ad114-7e13-5084-954f-fe012c677804"
 
+[[MozillaCACerts_jll]]
+uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
+version = "2022.10.11"
+
+[[NetworkOptions]]
+uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
+version = "1.2.0"
+
+[[OffsetArrays]]
+git-tree-sha1 = "1a27764e945a152f7ca7efa04de513d473e9542e"
+uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
+version = "1.14.1"
+
+    [OffsetArrays.extensions]
+    OffsetArraysAdaptExt = "Adapt"
+
+    [OffsetArrays.weakdeps]
+    Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
+
+[[OpenBLAS_jll]]
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
+uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
+version = "0.3.21+4"
+
 [[OrderedCollections]]
 deps = ["Random", "Serialization", "Test"]
 git-tree-sha1 = "c4c13474d23c60d20a67b217f1d7f22a40edf8f1"
@@ -117,8 +190,9 @@ uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
 version = "0.3.10"
 
 [[Pkg]]
-deps = ["Dates", "LibGit2", "Markdown", "Printf", "REPL", "Random", "SHA", "UUIDs"]
+deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
 uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+version = "1.9.2"
 
 [[Printf]]
 deps = ["Unicode"]
@@ -137,11 +211,11 @@ uuid = "d330b81b-6aea-500a-939a-2ce795aea3ee"
 version = "2.8.2"
 
 [[REPL]]
-deps = ["InteractiveUtils", "Markdown", "Sockets"]
+deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
 uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
 
 [[Random]]
-deps = ["Serialization"]
+deps = ["SHA", "Serialization"]
 uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
 [[Reexport]]
@@ -152,6 +226,7 @@ version = "0.2.0"
 
 [[SHA]]
 uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
+version = "0.7.0"
 
 [[Serialization]]
 uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
@@ -164,15 +239,31 @@ uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"
 uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
 
 [[SparseArrays]]
-deps = ["LinearAlgebra", "Random"]
+deps = ["Libdl", "LinearAlgebra", "Random", "Serialization", "SuiteSparse_jll"]
 uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [[Statistics]]
 deps = ["LinearAlgebra", "SparseArrays"]
 uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+version = "1.9.0"
+
+[[SuiteSparse_jll]]
+deps = ["Artifacts", "Libdl", "Pkg", "libblastrampoline_jll"]
+uuid = "bea87d4a-7f5b-5778-9afe-8cc45184846c"
+version = "5.10.1+6"
+
+[[TOML]]
+deps = ["Dates"]
+uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
+version = "1.0.3"
+
+[[Tar]]
+deps = ["ArgTools", "SHA"]
+uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"
+version = "1.10.0"
 
 [[Test]]
-deps = ["Distributed", "InteractiveUtils", "Logging", "Random"]
+deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
 uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [[UUIDs]]
@@ -186,3 +277,23 @@ uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
 git-tree-sha1 = "80229be1f670524750d905f8fc8148e5a8c4537f"
 uuid = "81def892-9a0e-5fdd-b105-ffc91e053289"
 version = "1.2.0"
+
+[[Zlib_jll]]
+deps = ["Libdl"]
+uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
+version = "1.2.13+0"
+
+[[libblastrampoline_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
+version = "5.8.0+0"
+
+[[nghttp2_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"
+version = "1.52.0+1"
+
+[[p7zip_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
+version = "17.4.0+0"

docs/Project.toml

@@ -1,5 +1,6 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+FLOWMath = "6cb5d3fb-0fe8-4cc2-bd89-9fe0b19a99d3"
 PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
 
 [compat]

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

@@ -15,4 +15,36 @@ function trapz(x, y)
         integral += (x[i+1]-x[i])*0.5*(y[i] + y[i+1])
     end
     return integral
-end
+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`.
+"""
+function cumtrapz(x, y)
+    integral = similar(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
+
+
+"""
+    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 ------