nathanaelbosch / ProbNumDiffEq.jl

Probabilistic Numerical Differential Equation solvers via Bayesian filtering and smoothing
MIT License
119 stars 16 forks source link

Use Gaussian as Initial Condition #318

Closed TheFibonacciEffect closed 3 months ago

TheFibonacciEffect commented 3 months ago

Is it possible to use a Gaussian as an intial condition?

For example something like this?

using LinearAlgebra, Statistics, Distributions
using DiffEqDevTools, ParameterizedFunctions, SciMLBase, OrdinaryDiffEq, Sundials, Plots, ODEInterfaceDiffEq
using ModelingToolkit
using ProbNumDiffEq

function osscilator!(ddu, du, u, p,t)
    ddu .= - p .* u
end

# osscilator
u0 = [1]
du0 = [1]
p = 100
T = 3
t = 0:0.1:T
prob = SecondOrderODEProblem(osscilator!, du0, u0, (0,T),p)
@time sol = solve(prob, EK0(;smooth=true), abstol=1e-1, reltol=1e-1)
plot(sol)

Where I would like to use u0 = [Gaussian(1,1)] or something similar as an intial condition instead.

Currently it does not seem to be suported: ERROR: MethodError: no method matching zero(::Type{Any}) but it would be very useful. For example when the initial condition is a result of a measurement and is not known with infinite precision.

nathanaelbosch commented 3 months ago

Thanks for raising this issue! This is a common question I encounter. Unfortunately probabilistic ODE solvers cannot actually handle Gaussian initial conditions, which is why this functionality is also not implemented.

To give a short explanation, essentially filters and smoothers do not "propagate" uncertainty but rather "estimate" the unknown states. By this I mean that if we start with an uncertain initial value, the smoother ends up computing a posterior distribution over this initial value, such that it leads to a higher likelihood of satisfying the ODE. But this is most likely not what you're after. Instead, the problem you are describing is one where you want to keep exactly the uncertainty that you specified initially, and propagate it to the resulting state trajectory (by marginalizing it out). There is currently no published specialied approach for how to do this within the framework of filtering-based probabilistic ODE solvers, so your best bet is to use any standard uncetainty propagation method (e.g. sampling) on top of a probabilistic numerical method.

I hope this clarifies your question! I will close this issue for now, but if you have any follow-up question I'm also happy to continue the discussion.