From 45d7aa39ccecc77684169db9e6f5e64a5226a6e4 Mon Sep 17 00:00:00 2001 From: Daniel Ingraham Date: Fri, 16 Aug 2024 11:39:34 -0400 Subject: [PATCH] Add first try at `cumtrapz`/`cumtrapz!` --- docs/Manifest.toml | 123 ++++++++++++++++++++++++++++++++++++++++++--- docs/Project.toml | 1 + docs/src/index.md | 25 +++++++++ src/FLOWMath.jl | 2 +- src/quadrature.jl | 34 ++++++++++++- test/runtests.jl | 36 ++++++++++++- 6 files changed, 212 insertions(+), 9 deletions(-) diff --git a/docs/Manifest.toml b/docs/Manifest.toml index a11b446..fdba91f 100644 --- a/docs/Manifest.toml +++ b/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" diff --git a/docs/Project.toml b/docs/Project.toml index 7c05cc0..e47fd2d 100644 --- a/docs/Project.toml +++ b/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] diff --git a/docs/src/index.md b/docs/src/index.md index fee7c47..df0c7a3 100644 --- a/docs/src/index.md +++ b/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) diff --git a/src/FLOWMath.jl b/src/FLOWMath.jl index b5006aa..3bcb5e4 100644 --- a/src/FLOWMath.jl +++ b/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 diff --git a/src/quadrature.jl b/src/quadrature.jl index 022fd79..b1f8d14 100644 --- a/src/quadrature.jl +++ b/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 \ No newline at end of file +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 diff --git a/test/runtests.jl b/test/runtests.jl index 1c7d39c..44fd711 100644 --- a/test/runtests.jl +++ b/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 ------