The problem is somehow method-dependent. E.g. when using Trapezoid() the code
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
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