# Why does my cubature of ODESolution break down with finite limits?

**URL:** <https://discourse.julialang.org/t/why-does-my-cubature-of-odesolution-break-down-with-finite-limits/108736>\
**Category:** Numerics\
**Tags:** diffeq, distributions, integral, differentialequation, ordinarydiffeq\
**Created:** [January 12, 2024, 4:25pm UTC](https://discourse.julialang.org/t/why-does-my-cubature-of-odesolution-break-down-with-finite-limits/108736 "2024-01-12T16:25:36Z")\
**Posts on this page:** 1\
**Showing post:** 8

<div class="post-metadata">

**Author:** ![homocomputeris](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/homocomputeris/32/8933_2.png) [@homocomputeris](https://discourse.julialang.org/u/homocomputeris)\
**Post date:** [January 16, 2024, 9:47pm UTC](https://discourse.julialang.org/t/why-does-my-cubature-of-odesolution-break-down-with-finite-limits/108736/8 "2024-01-16T21:47:15Z")

</div>

The problem is somehow method-dependent. E.g. when using `Trapezoid()` the code

```julia
module Mwe

using DifferentialEquations
using Integrals
using Distributions
using LinearAlgebra
using ForwardDiff

# just to speed things up a bit
const TOL = 1e-1

# some parameters; exact values do not matter
tspan = (0, pi + 0.0)
p = [1.0, 1.0]
t = 22.0 / 7.0
# limit the domain of interest
rectangle = (-5, 5)

function f(u, p, t)
    du = [u[2], u[1], 0.0]
    return du
end

# scalar function with vector argument
function g0(x, p, t)
    return pdf(truncated(Normal(1, 1), rectangle...), x[1]) *
           pdf(truncated(Normal(1, 1), rectangle...), x[2]) *
           pdf(truncated(Normal(1, 1), rectangle...), x[3])
end

function g1(x, p, t)
    # find the initial condition for a given point
    function tmp(y)
        prob = ODEProblem(f, y, reverse(tspan), p)
        sol = solve(prob, Trapezoid(), dt = pi / 10)
        return sol(tspan[end] - t)
    end
    # use the found IC in g0
    return g0(tmp(x), p, 0.0) * abs(det(ForwardDiff.jacobian(tmp, x)))
end

# Truncate in a rectangle
function x1_broken(jpdf, x::Number, p, t)
    prob = IntegralProblem((tail, p) -> jpdf([x; tail], p, t),
        [-5, rectangle[1]],
        [5, rectangle[2]], p)
    sol = solve(prob, HCubatureJL();
        reltol = TOL, abstol = TOL)
    return sol.u
end

@show g0(rand(3), p, t)
@show g1(rand(3), p, t)

# this doesn't
@show x1_broken(g1, 3.0, p, t)

end

```

outputs

```julia
g0(rand(3), p, t) = 0.03282973041989611
g1(rand(3), p, t) = 0.004435854935625869
ERROR: MethodError: Cannot `convert` an object of type LinearAlgebra.LU{Float64, Matrix{Float64}, Vector{Int64}} to an object of type StaticArrays.LU{LinearAlgebra.LowerTriangular{Float64, StaticArraysCore.SMatrix{3, 3, Float64, 9}}, LinearAlgebra.UpperTriangular{Float64, StaticArraysCore.SMatrix{3, 3, Float64, 9}}, StaticArraysCore.SVector{3, Int64}}

Closest candidates are:
  convert(::Type{T}, ::T) where T
   @ Base Base.jl:84
  (::Type{StaticArrays.LU{L, U, p}} where {L, U, p})(::Any, ::Any, ::Any)
   @ StaticArrays ~/.julia/packages/StaticArrays/PLKkM/src/lu.jl:3

```

---

_[View the full topic](https://discourse.julialang.org/t/why-does-my-cubature-of-odesolution-break-down-with-finite-limits/108736)._