diff --git a/docs/make.jl b/docs/make.jl
index 8a2835d7..94546761 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -13,6 +13,6 @@ makedocs(modules = [DataInterpolations],
     format = Documenter.HTML(assets = ["assets/favicon.ico"],
         canonical = "https://docs.sciml.ai/DataInterpolations/stable/"),
     pages = ["index.md", "Methods" => "methods.md",
-        "Interface" => "interface.md", "Manual" => "manual.md"])
+        "Interface" => "interface.md", "Manual" => "manual.md", "Inverting Integrals" => "inverting_integrals.md"])
 
 deploydocs(repo = "github.com/SciML/DataInterpolations.jl"; push_preview = true)
diff --git a/docs/src/inverting_integrals.md b/docs/src/inverting_integrals.md
new file mode 100644
index 00000000..8cbfcc3e
--- /dev/null
+++ b/docs/src/inverting_integrals.md
@@ -0,0 +1,52 @@
+# Inverting integrals
+
+Solving implicit integral problems of the form:
+
+```math
+\begin{equation}
+    \text{find $t$ such that } \int_{t_1}^t f(\tau)\text{d}\tau = V \ge 0
+\end{equation}
+```
+
+is supported for interpolations $f$ that are strictly positive and of one of these types:
+
+  - `ConstantInterpolation`
+  - `LinearInterpolation`
+
+This is done by creating an 'integral inverse' interpolation object which can efficiently compute $t$ for a given value of $V$, see the example below.
+
+```@example inverting_integrals
+using Random #hide
+Random.seed!(1234) # hide
+using DataInterpolations
+using Plots
+
+# Create LinearInterpolation object from the
+u = sqrt.(1:25) + (2.0 * rand(25) .- 1.0) / 3
+t = cumsum(rand(25))
+A = LinearInterpolation(u, t)
+
+# Create LinearInterpolationIntInv object
+# from the LinearInterpolation object
+A_intinv = DataInterpolations.invert_integral(A)
+
+# Get the t values up to and including the
+# solution to the integral problem
+V = 25.0
+t_ = A_intinv(V)
+ts = t[t .<= t_]
+push!(ts, t_)
+
+# Plot results
+plot(A; label = "Linear Interpolation")
+plot!(ts, A.(ts), fillrange = 0.0, fillalpha = 0.75,
+    fc = :blues, lw = 0, label = "Area of $V")
+```
+
+## Docstrings
+
+```@docs
+DataInterpolations.invert_integral
+ConstantInterpolationIntInv
+LinearInterpolationIntInv
+```
diff --git a/src/DataInterpolations.jl b/src/DataInterpolations.jl
index 9df9db68..c86a6579 100644
--- a/src/DataInterpolations.jl
+++ b/src/DataInterpolations.jl
@@ -18,6 +18,7 @@ include("interpolation_methods.jl")
 include("plot_rec.jl")
 include("derivatives.jl")
 include("integrals.jl")
+include("integral_inverses.jl")
 include("online.jl")
 include("show.jl")
 
@@ -75,6 +76,18 @@ function Base.showerror(io::IO, e::DerivativeNotFoundError)
     print(io, DERIVATIVE_NOT_FOUND_ERROR)
 end
 
+const INTEGRAL_INVERSE_NOT_FOUND_ERROR = "Cannot invert the integral analytically. Please use Numerical methods."
+struct IntegralInverseNotFoundError <: Exception end
+function Base.showerror(io::IO, e::IntegralInverseNotFoundError)
+    print(io, INTEGRAL_INVERSE_NOT_FOUND_ERROR)
+end
+
+const INTEGRAL_NOT_INVERTIBLE_ERROR = "The Interpolation is not positive everywhere so its integral is not invertible."
+struct IntegralNotInvertibleError <: Exception end
+function Base.showerror(io::IO, e::IntegralNotInvertibleError)
+    print(io, INTEGRAL_NOT_INVERTIBLE_ERROR)
+end
+
 const MUST_COPY_ERROR = "A copy must be made of u, t to filter missing data"
 struct MustCopyError <: Exception end
 function Base.showerror(io::IO, e::MustCopyError)
@@ -84,7 +97,7 @@ end
 export LinearInterpolation, QuadraticInterpolation, LagrangeInterpolation,
        AkimaInterpolation, ConstantInterpolation, QuadraticSpline, CubicSpline,
        BSplineInterpolation, BSplineApprox, CubicHermiteSpline,
-       QuinticHermiteSpline
+       QuinticHermiteSpline, LinearInterpolationIntInv, ConstantInterpolationIntInv
 
 # added for RegularizationSmooth, JJS 11/27/21
 ### Regularization data smoothing and interpolation
diff --git a/src/integral_inverses.jl b/src/integral_inverses.jl
new file mode 100644
index 00000000..1e2ca581
--- /dev/null
+++ b/src/integral_inverses.jl
@@ -0,0 +1,95 @@
+abstract type AbstractIntegralInverseInterpolation{T} <: AbstractInterpolation{T} end
+
+"""
+    invert_integral(A::AbstractInterpolation)::AbstractIntegralInverseInterpolation
+
+Creates the inverted integral interpolation object from the given interpolation. Conditions:
+
+  - The range of `A` must be strictly positive
+  - There must be an ordering defined on the data type of `A.u`
+  - This is currently only supported for ConstantInterpolation and LinearInterpolation
+
+## Arguments
+
+  - `A`: interpolation object satisfying the above requirements
+"""
+invert_integral(A::AbstractInterpolation) = throw(IntegralInverseNotFoundError())
+
+_integral(A::AbstractIntegralInverseInterpolation, idx, t) = throw(IntegralNotFoundError())
+
+function _derivative(A::AbstractIntegralInverseInterpolation, t::Number, iguess)
+    inv(A.itp(A(t))), A.idx_prev[]
+end
+
+"""
+    some stuff
+"""
+struct LinearInterpolationIntInv{uType, tType, itpType, T} <:
+       AbstractIntegralInverseInterpolation{T}
+    u::uType
+    t::tType
+    extrapolate::Bool
+    idx_prev::Base.RefValue{Int}
+    itp::itpType
+    function LinearInterpolationIntInv(u, t, A)
+        new{typeof(u), typeof(t), typeof(A), eltype(u)}(
+            u, t, A.extrapolate, Ref(1), A)
+    end
+end
+
+function invertible_integral(A::LinearInterpolation{<:AbstractVector{<:Number}})
+    return all(A.u .> 0)
+end
+
+function invert_integral(A::LinearInterpolation{<:AbstractVector{<:Number}})
+    !invertible_integral(A) && throw(IntegralNotInvertibleError())
+    return LinearInterpolationIntInv(A.t, A.I, A)
+end
+
+function _interpolate(
+        A::LinearInterpolationIntInv{<:AbstractVector{<:Number}}, t::Number, iguess)
+    idx = get_idx(A.t, t, iguess)
+    Δt = t - A.t[idx]
+    x = A.itp.u[idx]
+    u = A.u[idx] + 2Δt / (x + sqrt(x^2 + A.itp.p.slope[idx] * 2Δt))
+    u, idx
+end
+
+"""
+    some stuff
+"""
+struct ConstantInterpolationIntInv{uType, tType, itpType, T} <:
+       AbstractIntegralInverseInterpolation{T}
+    u::uType
+    t::tType
+    extrapolate::Bool
+    idx_prev::Base.RefValue{Int}
+    itp::itpType
+    function ConstantInterpolationIntInv(u, t, A)
+        new{typeof(u), typeof(t), typeof(A), eltype(u)}(
+            u, t, A.extrapolate, Ref(1), A
+        )
+    end
+end
+
+function invertible_integral(A::ConstantInterpolation{<:AbstractVector{<:Number}})
+    return all(A.u .> 0)
+end
+
+function invert_integral(A::ConstantInterpolation{<:AbstractVector{<:Number}})
+    !invertible_integral(A) && throw(IntegralNotInvertibleError())
+    return ConstantInterpolationIntInv(A.t, A.I, A)
+end
+
+function _interpolate(
+        A::ConstantInterpolationIntInv{<:AbstractVector{<:Number}}, t::Number, iguess)
+    idx = get_idx(A.t, t, iguess; ub_shift = 0)
+    if A.itp.dir === :left
+        # :left means that value to the left is used for interpolation
+        idx_ = get_idx(A.t, t, idx; lb = 1, ub_shift = 0)
+    else
+        # :right means that value to the right is used for interpolation
+        idx_ = get_idx(A.t, t, idx; side = :first, lb = 1, ub_shift = 0)
+    end
+    A.u[idx] + (t - A.t[idx]) / A.itp.u[idx_], idx
+end
diff --git a/src/integrals.jl b/src/integrals.jl
index e64b7568..3040189f 100644
--- a/src/integrals.jl
+++ b/src/integrals.jl
@@ -6,20 +6,22 @@ end
 function integral(A::AbstractInterpolation, t1::Number, t2::Number)
     ((t1 < A.t[1] || t1 > A.t[end]) && !A.extrapolate) && throw(ExtrapolationError())
     ((t2 < A.t[1] || t2 > A.t[end]) && !A.extrapolate) && throw(ExtrapolationError())
+    !hasfield(typeof(A), :I) && throw(IntegralNotFoundError())
     # the index less than or equal to t1
-    idx1 = max(1, min(searchsortedlast(A.t, t1), length(A.t) - 1))
+    idx1 = get_idx(A.t, t1, 0)
     # the index less than t2
-    idx2 = max(1, min(searchsortedlast(A.t, t2), length(A.t) - 1))
-    if A.t[idx2] == t2
-        idx2 -= 1
-    end
-    total = zero(eltype(A.u))
-    for idx in idx1:idx2
-        lt1 = idx == idx1 ? t1 : A.t[idx]
-        lt2 = idx == idx2 ? t2 : A.t[idx + 1]
-        total += _integral(A, idx, lt2) - _integral(A, idx, lt1)
+    idx2 = get_idx(A.t, t2, 0; idx_shift = -1, side = :first)
+
+    total = A.I[idx2] - A.I[idx1]
+    return if t1 == t2
+        zero(total)
+    else
+        total += _integral(A, idx1, A.t[idx1])
+        total -= _integral(A, idx1, t1)
+        total += _integral(A, idx2, t2)
+        total -= _integral(A, idx2, A.t[idx2])
+        total
     end
-    total
 end
 
 function _integral(A::LinearInterpolation{<:AbstractVector{<:Number}},
@@ -29,7 +31,8 @@ function _integral(A::LinearInterpolation{<:AbstractVector{<:Number}},
     Δt * (A.u[idx] + A.p.slope[idx] * Δt / 2)
 end
 
-function _integral(A::ConstantInterpolation{<:AbstractVector}, idx::Number, t::Number)
+function _integral(
+        A::ConstantInterpolation{<:AbstractVector{<:Number}}, idx::Number, t::Number)
     if A.dir === :left
         # :left means that value to the left is used for interpolation
         return A.u[idx] * t
diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl
index bb6339fe..83e04fe5 100644
--- a/src/interpolation_caches.jl
+++ b/src/interpolation_caches.jl
@@ -14,23 +14,26 @@ Extrapolation extends the last linear polynomial on each side.
   - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
   - `safetycopy`: boolean value to make a copy of `u` and `t`. Defaults to `true`.
 """
-struct LinearInterpolation{uType, tType, pType, T} <: AbstractInterpolation{T}
+struct LinearInterpolation{uType, tType, IType, pType, T} <: AbstractInterpolation{T}
     u::uType
     t::tType
+    I::IType
     p::LinearParameterCache{pType}
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
     safetycopy::Bool
-    function LinearInterpolation(u, t, p, extrapolate, safetycopy)
-        new{typeof(u), typeof(t), typeof(p.slope), eltype(u)}(
-            u, t, p, extrapolate, Ref(1), safetycopy)
+    function LinearInterpolation(u, t, I, p, extrapolate, safetycopy)
+        new{typeof(u), typeof(t), typeof(I), typeof(p.slope), eltype(u)}(
+            u, t, I, p, extrapolate, Ref(1), safetycopy)
     end
 end
 
 function LinearInterpolation(u, t; extrapolate = false, safetycopy = true)
     u, t = munge_data(u, t, safetycopy)
     p = LinearParameterCache(u, t)
-    LinearInterpolation(u, t, p, extrapolate, safetycopy)
+    A = LinearInterpolation(u, t, nothing, p, extrapolate, safetycopy)
+    I = cumulative_integral(A)
+    LinearInterpolation(u, t, I, p, extrapolate, safetycopy)
 end
 
 """
@@ -50,32 +53,33 @@ Extrapolation extends the last quadratic polynomial on each side.
   - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
   - `safetycopy`: boolean value to make a copy of `u` and `t`. Defaults to `true`.
 """
-struct QuadraticInterpolation{uType, tType, pType, T} <: AbstractInterpolation{T}
+struct QuadraticInterpolation{uType, tType, IType, pType, T} <: AbstractInterpolation{T}
     u::uType
     t::tType
+    I::IType
     p::QuadraticParameterCache{pType}
     mode::Symbol
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
     safetycopy::Bool
-    function QuadraticInterpolation(u, t, p, mode, extrapolate, safetycopy)
+    function QuadraticInterpolation(u, t, I, p, mode, extrapolate, safetycopy)
         mode ∈ (:Forward, :Backward) ||
             error("mode should be :Forward or :Backward for QuadraticInterpolation")
-        new{typeof(u), typeof(t), typeof(p.l₀), eltype(u)}(
-            u, t, p, mode, extrapolate, Ref(1), safetycopy)
+        new{typeof(u), typeof(t), typeof(I), typeof(p.l₀), eltype(u)}(
+            u, t, I, p, mode, extrapolate, Ref(1), safetycopy)
     end
 end
 
 function QuadraticInterpolation(u, t, mode; extrapolate = false, safetycopy = true)
     u, t = munge_data(u, t, safetycopy)
     p = QuadraticParameterCache(u, t)
-    QuadraticInterpolation(u, t, p, mode, extrapolate, safetycopy)
+    A = QuadraticInterpolation(u, t, nothing, p, mode, extrapolate, safetycopy)
+    I = cumulative_integral(A)
+    QuadraticInterpolation(u, t, I, p, mode, extrapolate, safetycopy)
 end
 
 function QuadraticInterpolation(u, t; extrapolate = false, safetycopy = true)
-    u, t = munge_data(u, t, safetycopy)
-    p = QuadraticParameterCache(u, t)
-    QuadraticInterpolation(u, t, p, :Forward, extrapolate, safetycopy)
+    QuadraticInterpolation(u, t, :Forward; extrapolate, safetycopy)
 end
 
 """
@@ -145,20 +149,22 @@ Extrapolation extends the last cubic polynomial on each side.
   - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
   - `safetycopy`: boolean value to make a copy of `u` and `t`. Defaults to `true`.
 """
-struct AkimaInterpolation{uType, tType, bType, cType, dType, T} <:
+struct AkimaInterpolation{uType, tType, IType, bType, cType, dType, T} <:
        AbstractInterpolation{T}
     u::uType
     t::tType
+    I::IType
     b::bType
     c::cType
     d::dType
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
     safetycopy::Bool
-    function AkimaInterpolation(u, t, b, c, d, extrapolate, safetycopy)
-        new{typeof(u), typeof(t), typeof(b), typeof(c),
+    function AkimaInterpolation(u, t, I, b, c, d, extrapolate, safetycopy)
+        new{typeof(u), typeof(t), typeof(I), typeof(b), typeof(c),
             typeof(d), eltype(u)}(u,
             t,
+            I,
             b,
             c,
             d,
@@ -191,7 +197,9 @@ function AkimaInterpolation(u, t; extrapolate = false, safetycopy = true)
     c = (3.0 .* m[3:(end - 2)] .- 2.0 .* b[1:(end - 1)] .- b[2:end]) ./ dt
     d = (b[1:(end - 1)] .+ b[2:end] .- 2.0 .* m[3:(end - 2)]) ./ dt .^ 2
 
-    AkimaInterpolation(u, t, b, c, d, extrapolate, safetycopy)
+    A = AkimaInterpolation(u, t, nothing, b, c, d, extrapolate, safetycopy)
+    I = cumulative_integral(A)
+    AkimaInterpolation(u, t, I, b, c, d, extrapolate, safetycopy)
 end
 
 """
@@ -212,23 +220,26 @@ Extrapolation extends the last constant polynomial at the end points on each sid
   - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
   - `safetycopy`: boolean value to make a copy of `u` and `t`. Defaults to `true`.
 """
-struct ConstantInterpolation{uType, tType, dirType, T} <: AbstractInterpolation{T}
+struct ConstantInterpolation{uType, tType, IType, T} <: AbstractInterpolation{T}
     u::uType
     t::tType
+    I::IType
     p::Nothing
     dir::Symbol # indicates if value to the $dir should be used for the interpolation
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
     safetycopy::Bool
-    function ConstantInterpolation(u, t, dir, extrapolate, safetycopy)
-        new{typeof(u), typeof(t), typeof(dir), eltype(u)}(
-            u, t, nothing, dir, extrapolate, Ref(1), safetycopy)
+    function ConstantInterpolation(u, t, I, dir, extrapolate, safetycopy)
+        new{typeof(u), typeof(t), typeof(I), eltype(u)}(
+            u, t, I, nothing, dir, extrapolate, Ref(1), safetycopy)
     end
 end
 
 function ConstantInterpolation(u, t; dir = :left, extrapolate = false, safetycopy = true)
     u, t = munge_data(u, t, safetycopy)
-    ConstantInterpolation(u, t, dir, extrapolate, safetycopy)
+    A = ConstantInterpolation(u, t, nothing, dir, extrapolate, safetycopy)
+    I = cumulative_integral(A)
+    ConstantInterpolation(u, t, I, dir, extrapolate, safetycopy)
 end
 
 """
@@ -247,10 +258,11 @@ Extrapolation extends the last quadratic polynomial on each side.
   - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
   - `safetycopy`: boolean value to make a copy of `u` and `t`. Defaults to `true`.
 """
-struct QuadraticSpline{uType, tType, pType, tAType, dType, zType, T} <:
+struct QuadraticSpline{uType, tType, IType, pType, tAType, dType, zType, T} <:
        AbstractInterpolation{T}
     u::uType
     t::tType
+    I::IType
     p::QuadraticSplineParameterCache{pType}
     tA::tAType
     d::dType
@@ -258,10 +270,11 @@ struct QuadraticSpline{uType, tType, pType, tAType, dType, zType, T} <:
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
     safetycopy::Bool
-    function QuadraticSpline(u, t, p, tA, d, z, extrapolate, safetycopy)
-        new{typeof(u), typeof(t), typeof(p.σ), typeof(tA),
+    function QuadraticSpline(u, t, I, p, tA, d, z, extrapolate, safetycopy)
+        new{typeof(u), typeof(t), typeof(I), typeof(p.σ), typeof(tA),
             typeof(d), typeof(z), eltype(u)}(u,
             t,
+            I,
             p,
             tA,
             d,
@@ -273,9 +286,9 @@ struct QuadraticSpline{uType, tType, pType, tAType, dType, zType, T} <:
     end
 end
 
-function QuadraticSpline(u::uType,
-        t;
-        extrapolate = false, safetycopy = true) where {uType <: AbstractVector{<:Number}}
+function QuadraticSpline(
+        u::uType, t; extrapolate = false,
+        safetycopy = true) where {uType <: AbstractVector{<:Number}}
     u, t = munge_data(u, t, safetycopy)
     s = length(t)
     dl = ones(eltype(t), s - 1)
@@ -289,7 +302,9 @@ function QuadraticSpline(u::uType,
     d = map(i -> i == 1 ? typed_zero : 2 // 1 * (u[i] - u[i - 1]) / (t[i] - t[i - 1]), 1:s)
     z = tA \ d
     p = QuadraticSplineParameterCache(z, t)
-    QuadraticSpline(u, t, p, tA, d, z, extrapolate, safetycopy)
+    A = QuadraticSpline(u, t, nothing, p, tA, d, z, extrapolate, safetycopy)
+    I = cumulative_integral(A)
+    QuadraticSpline(u, t, I, p, tA, d, z, extrapolate, safetycopy)
 end
 
 function QuadraticSpline(
@@ -308,7 +323,9 @@ function QuadraticSpline(
     z_ = reshape(transpose(tA \ d), size(u[1])..., :)
     z = [z_s for z_s in eachslice(z_, dims = ndims(z_))]
     p = QuadraticSplineParameterCache(z, t)
-    QuadraticSpline(u, t, p, tA, d, z, extrapolate, safetycopy)
+    A = QuadraticSpline(u, t, nothing, p, tA, d, z, extrapolate, safetycopy)
+    I = cumulative_integral(A)
+    QuadraticSpline(u, t, I, p, tA, d, z, extrapolate, safetycopy)
 end
 
 """
@@ -327,18 +344,21 @@ Second derivative on both ends are zero, which are also called "natural" boundar
   - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
   - `safetycopy`: boolean value to make a copy of `u` and `t`. Defaults to `true`.
 """
-struct CubicSpline{uType, tType, pType, hType, zType, T} <: AbstractInterpolation{T}
+struct CubicSpline{uType, tType, IType, pType, hType, zType, T} <: AbstractInterpolation{T}
     u::uType
     t::tType
+    I::IType
     p::CubicSplineParameterCache{pType}
     h::hType
     z::zType
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
     safetycopy::Bool
-    function CubicSpline(u, t, p, h, z, extrapolate, safetycopy)
-        new{typeof(u), typeof(t), typeof(p.c₁), typeof(h), typeof(z), eltype(u)}(u,
+    function CubicSpline(u, t, I, p, h, z, extrapolate, safetycopy)
+        new{typeof(u), typeof(t), typeof(I), typeof(p.c₁), typeof(h), typeof(z), eltype(u)}(
+            u,
             t,
+            I,
             p,
             h,
             z,
@@ -370,7 +390,9 @@ function CubicSpline(u::uType,
         1:(n + 1))
     z = tA \ d
     p = CubicSplineParameterCache(u, h, z)
-    CubicSpline(u, t, p, h[1:(n + 1)], z, extrapolate, safetycopy)
+    A = CubicSpline(u, t, nothing, p, h[1:(n + 1)], z, extrapolate, safetycopy)
+    I = cumulative_integral(A)
+    CubicSpline(u, t, I, p, h[1:(n + 1)], z, extrapolate, safetycopy)
 end
 
 function CubicSpline(
@@ -390,7 +412,9 @@ function CubicSpline(
     z_ = reshape(transpose(tA \ d), size(u[1])..., :)
     z = [z_s for z_s in eachslice(z_, dims = ndims(z_))]
     p = CubicSplineParameterCache(u, h, z)
-    CubicSpline(u, t, p, h[1:(n + 1)], z, extrapolate, safetycopy)
+    A = CubicSpline(u, t, nothing, p, h[1:(n + 1)], z, extrapolate, safetycopy)
+    I = cumulative_integral(A)
+    CubicSpline(u, t, I, p, h[1:(n + 1)], z, extrapolate, safetycopy)
 end
 
 """
@@ -691,17 +715,18 @@ It is a Cubic Hermite interpolation, which is a piece-wise third degree polynomi
   - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
   - `safetycopy`: boolean value to make a copy of `u` and `t`. Defaults to `true`.
 """
-struct CubicHermiteSpline{uType, tType, duType, pType, T} <: AbstractInterpolation{T}
+struct CubicHermiteSpline{uType, tType, IType, duType, pType, T} <: AbstractInterpolation{T}
     du::duType
     u::uType
     t::tType
+    I::IType
     p::CubicHermiteParameterCache{pType}
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
     safetycopy::Bool
-    function CubicHermiteSpline(du, u, t, p, extrapolate, safetycopy)
-        new{typeof(u), typeof(t), typeof(du), typeof(p.c₁), eltype(u)}(
-            du, u, t, p, extrapolate, Ref(1), safetycopy)
+    function CubicHermiteSpline(du, u, t, I, p, extrapolate, safetycopy)
+        new{typeof(u), typeof(t), typeof(I), typeof(du), typeof(p.c₁), eltype(u)}(
+            du, u, t, I, p, extrapolate, Ref(1), safetycopy)
     end
 end
 
@@ -709,7 +734,9 @@ function CubicHermiteSpline(du, u, t; extrapolate = false, safetycopy = true)
     @assert length(u)==length(du) "Length of `u` is not equal to length of `du`."
     u, t = munge_data(u, t, safetycopy)
     p = CubicHermiteParameterCache(du, u, t)
-    return CubicHermiteSpline(du, u, t, p, extrapolate, safetycopy)
+    A = CubicHermiteSpline(du, u, t, nothing, p, extrapolate, safetycopy)
+    I = cumulative_integral(A)
+    CubicHermiteSpline(du, u, t, I, p, extrapolate, safetycopy)
 end
 
 """
@@ -729,19 +756,21 @@ It is a Quintic Hermite interpolation, which is a piece-wise fifth degree polyno
   - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
   - `safetycopy`: boolean value to make a copy of `u` and `t`. Defaults to `true`.
 """
-struct QuinticHermiteSpline{uType, tType, duType, dduType, pType, T} <:
+struct QuinticHermiteSpline{uType, tType, IType, duType, dduType, pType, T} <:
        AbstractInterpolation{T}
     ddu::dduType
     du::duType
     u::uType
     t::tType
+    I::IType
     p::QuinticHermiteParameterCache{pType}
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
     safetycopy::Bool
-    function QuinticHermiteSpline(ddu, du, u, t, p, extrapolate, safetycopy)
-        new{typeof(u), typeof(t), typeof(du), typeof(ddu), typeof(p.c₁), eltype(u)}(
-            ddu, du, u, t, p, extrapolate, Ref(1), safetycopy)
+    function QuinticHermiteSpline(ddu, du, u, t, I, p, extrapolate, safetycopy)
+        new{typeof(u), typeof(t), typeof(I), typeof(du),
+            typeof(ddu), typeof(p.c₁), eltype(u)}(
+            ddu, du, u, t, I, p, extrapolate, Ref(1), safetycopy)
     end
 end
 
@@ -749,5 +778,7 @@ function QuinticHermiteSpline(ddu, du, u, t; extrapolate = false, safetycopy = t
     @assert length(u)==length(du)==length(ddu) "Length of `u` is not equal to length of `du` or `ddu`."
     u, t = munge_data(u, t, safetycopy)
     p = QuinticHermiteParameterCache(ddu, du, u, t)
-    return QuinticHermiteSpline(ddu, du, u, t, p, extrapolate, safetycopy)
+    A = QuinticHermiteSpline(ddu, du, u, t, nothing, p, extrapolate, safetycopy)
+    I = cumulative_integral(A)
+    QuinticHermiteSpline(ddu, du, u, t, I, p, extrapolate, safetycopy)
 end
diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl
index 370c4062..466248b1 100644
--- a/src/interpolation_utils.jl
+++ b/src/interpolation_utils.jl
@@ -120,3 +120,13 @@ function get_idx(tvec, t, iguess; lb = 1, ub_shift = -1, idx_shift = 0, side = :
         error("side must be :first or :last")
     end
 end
+
+function cumulative_integral(A)
+    if isempty(methods(_integral, (typeof(A), Any, Any)))
+        return nothing
+    end
+    integral_values = [_integral(A, idx, A.t[idx + 1]) - _integral(A, idx, A.t[idx])
+                       for idx in 1:(length(A.t) - 1)]
+    pushfirst!(integral_values, zero(first(integral_values)))
+    return cumsum(integral_values)
+end
diff --git a/test/integral_inverse_tests.jl b/test/integral_inverse_tests.jl
new file mode 100644
index 00000000..fd83be67
--- /dev/null
+++ b/test/integral_inverse_tests.jl
@@ -0,0 +1,46 @@
+using DataInterpolations
+using DataInterpolations: integral, derivative, invert_integral
+using FiniteDifferences
+
+function test_integral_inverses(method; args = [], kwargs = [])
+    A = method(args...; kwargs..., extrapolate = true)
+    @test hasfield(typeof(A), :I)
+    A_intinv = invert_integral(A)
+    @test A_intinv isa DataInterpolations.AbstractIntegralInverseInterpolation
+    ts = range(first(A.t), last(A.t), length = 100)
+    Is = integral.(Ref(A), ts)
+    ts_ = A_intinv.(Is)
+    @test ts ≈ ts_
+
+    for I in Is
+        cdiff = forward_fdm(5, 1; geom = true)(A_intinv, I)
+        adiff = derivative(A_intinv, I)
+        @test cdiff ≈ adiff
+    end
+end
+
+@testset "Linear Interpolation" begin
+    t = collect(1:5)
+    u = [1.0, 1.0, 2.0, 4.0, 3.0]
+    test_integral_inverses(LinearInterpolation; args = [u, t])
+
+    u = [1.0, -1.0, 2.0, 4.0, 3.0]
+    A = LinearInterpolation(u, t)
+    @test_throws DataInterpolations.IntegralNotInvertibleError invert_integral(A)
+end
+
+@testset "Constant Interpolation" begin
+    t = collect(1:5)
+    u = [1.0, 1.0, 2.0, 4.0, 3.0]
+    test_integral_inverses(ConstantInterpolation; args = [u, t])
+    test_integral_inverses(ConstantInterpolation; args = [u, t], kwargs = [:dir => :right])
+
+    u = [1.0, -1.0, 2.0, 4.0, 3.0]
+    A = ConstantInterpolation(u, t)
+    @test_throws DataInterpolations.IntegralNotInvertibleError invert_integral(A)
+end
+
+t = collect(1:5)
+u = [1.0, 1.0, 2.0, 4.0, 3.0]
+A = QuadraticInterpolation(u, t)
+@test_throws DataInterpolations.IntegralInverseNotFoundError invert_integral(A)
diff --git a/test/runtests.jl b/test/runtests.jl
index 01f3a07e..0c722b2d 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -6,6 +6,7 @@ using SafeTestsets
 @safetestset "Interpolation Tests" include("interpolation_tests.jl")
 @safetestset "Derivative Tests" include("derivative_tests.jl")
 @safetestset "Integral Tests" include("integral_tests.jl")
+@safetestset "Integral Inverse Tests" include("integral_inverse_tests.jl")
 @safetestset "Online Tests" include("online_tests.jl")
 @safetestset "Regularization Smoothing" include("regularization.jl")
 @safetestset "Show methods" include("show.jl")