Validated Integration of Ordinary Differential Equations

In this last blog post, I'm going to talk about how to rigorously solve ordinary differential equations. These ideas are the essence of what is implemented in the PR#84. This is a work in progress that re-implement the validated ordinary differential equation solver and will continue even after the gsoc.

Ordinary Differential Equations (ODEs)

Let's define the problem, we look for a function that satisfies the following

Two remarks about the previous statement:

• We seek solutions not just for a one initial condition inline_formula not implemented but for a set (or box) of initial conditions.

• The solutions must be rigorous, this means that we can assert that the true solution of the ode is enclosed in some interval.

From the second point is clear that a good candidate to achieve this task are the Taylor Models.

First, note that every solution to the differential equation must satisfy the following equivalent integral equation

So, if we define the Picard (integral) operator as

we can restate the problem as the fixed point problem

The Schauder Theorem states the following:

Theorem: Let inline_formula not implemented be a continous operator on the Banach space inline_formula not implemented. Let inline_formula not implemented be compact and convex, and let inline_formula not implemented. Then inline_formula not implemented has a fixed point in inline_formula not implemented.

For our purposes, the operator inline_formula not implemented is the Picard integral operator and inline_formula not implemented is the space of continuous functions in inline_formula not implemented. Thus, to validate the solution all we have to do is to find a subset inline_formula not implemented such that inline_formula not implemented is entirely contained in inline_formula not implemented, and that the elements in inline_formula not implemented can be represented as Taylor Models.

Generalizing the Picard operator to Taylor Models, we get

This definition is the natural extension of the previous one, in this case, the integral is computed in a validated manner on Taylor Model arithmetic, and the initial Taylor Model inline_formula not implemented has the purpose of enclosing the set (or box) of initial conditions.

Therefore, to validate the solution of an ODE, we must find a Taylor Model such that

The critical point of the validated integration is: how do we actually find this Taylor Model?

Note that inline_formula not implemented is just the polynomial representation of the solution in a non-validated form, so, if we find inline_formula not implemented by some not-validated method (such as TaylorIntegration.jl), then we can just "fix it" by choosing an appropriate remainder inline_formula not implemented such that the last contention holds.

It turns out, that we can just take the polynomial part inline_formula not implementedand form an initial guess for the validated solution as inline_formula not implemented, then compute the Picard operator with this initial guess and iterate. In other words, iteratively we compute the sequence of Taylor Models

For some interation the fixed point equality must be fulfilled, and we can stop the iterative process.

The algorithm to compute a validated solution to the ODE problem can be summarized as:

1. Generate a Taylor Model such that the box of initial conditions at inline_formula not implemented is parametrized by it, this is inline_formula not implemented (for this time inline_formula not implemented the remainder inline_formula not implemented).

2. Given the time partition obtained from the ODE solver, do the following for each integration time step:

   1. Compute the polynomial part of the solution (by TaylorIntegration.jl), and iterate with the Picard operator until the contention inline_formula not implemented holds.

   2. Generate a new initial condition as inline_formula not implemented for the next integration step.

3. Repeat from 2.

To clarify this, let's go through a simple example.

]add TaylorModels IntervalArithmetic

using Plots

using TaylorModels

using Interact

using WebSockets

using WebIO

plot();

The example to look is the simple falling body problem, this is defined by the differential equation

let's take a box of initial conditions

# Defining the right hand side of the differential equation

@taylorize function fallingball!(dx, x, p, t)

    dx[1] = x[2]

    dx[2] = -p*one(x[1])

    nothing

end

g = 10. # gravity

symIbox = IntervalBox(-1 .. 1, Val(2))

q0 = [100.0, 0.0]

δq0 = IntervalBox(-5 .. 5, -0.5 .. 0.5)

tini, tend = 0.0, 4.

tstep = 1/16

tt = tini:tstep:tend

ξ = set_variables("ξₓ ξᵥ", order=2, numvars=length(q0))

# validated integration

tTM1, qv1, xTM1 = validated_integ(fallingball!, q0, δq0, tini, tend, 

    2, 4, 1.e-20, g, maxsteps=2_000);

The exact solution to this ode is

where inline_formula not implemented stands for the initial conditions.

# Exact solution

x(t, x0, g) = x0[1] + x0[2] * t - 0.5 * g * t^2

Let's plot the exact solution for a random selection of initial conditions and the validated solution for all solutions inside the initial condition box

interval_rand(X::Interval{T}) where {T} = X.lo + rand(T) * (X.hi - X.lo)

interval_rand(X::IntervalBox) = interval_rand.(X)

q0ξ = interval_rand(δq0)

print("Initial conditions x0 = $(q0ξ + q0)")

p = plot()

for i in 2:length(tTM1)

  t0 = tTM1[i-1]

  trange = range(0, xTM1[1, i].dom.hi, length=20)

  fT = xTM1[1, i]

  lo, hi = fT.rem.lo, fT.rem.hi

  mid_y = [mid(fT(t)(symIbox)) for t in trange]

  upper_y = ([fT(t)(symIbox).hi + hi for t in trange] - mid_y)

  lower_y = (mid_y - [fT(t)(symIbox).lo + abs(lo) for t in trange])

  plot!(trange.+t0, mid_y, ribbon=(lower_y, upper_y), 

    c=:dimgray, label="", lw=3)

end

plot!(tt, t -> x(t, q0ξ .+ q0, g), label="Exact solution", lw=2, c=:red2)

#plot!(tt, t -> x(t, q0 .+ mid.(δq0) .+ sup.(δq0), g), color=:green,

#  label="Sol inf")

#plot!(tt, t -> x(t, q0 .+ mid.(δq0) .+ inf.(δq0), g), color=:red,

#  label="Sol sup")

xlabel!("t")

ylabel!("x(t)")

p

Future work

Further improvements can be done to this validated ODE solver. In particular the step 2.1 from the previous algorithm can be improved by some additional procedure called inline_formula not implemented-inflation. This procedure explicitly modifies the remainder term as

in each application of the Picard operator. So, instead of just compute the remainder and iterate again, this remainder is modified in order to satisfy the inequality. This additional procedure can improve the convergence of the method to the validated Taylor Model.

This approach is going to be explored in the next couple of weeks, as can speed up the computation times of the current implementation.

Interactivity

If you want to play around and try different initial conditions you can use the following interactive plot.

tv = Dict()

tmv = Dict()

for δq01 = range(0, 10, length=10), δq02 = range(0, 10, length=10)

    q0 = [100.0, 0.0]

    δq0 = IntervalBox(-δq01 .. δq01, -δq02 .. δq02)

    tini, tend = 0.0, 4.

    tstep = 1/16

    tt = tini:tstep:tend

    ξ = set_variables("ξₓ ξᵥ", order=2, numvars=length(q0))

    # validated integration

    tTM1, qv1, xTM1 = validated_integ(fallingball!, q0, δq0, tini, tend, 

        2, 4, 1.e-20, g, maxsteps=2_000);

    tv["$(δq01)-$(δq02)"] = tTM1

    tmv["$(δq01)-$(δq02)"] = xTM1

end

@manipulate for δq01 = range(0, 10, length=10),

  δq02 = range(0, 10, length=10)

  q0 = [100.0, 0.0]

  δq0 = IntervalBox(-δq01 .. δq01, -δq02 .. δq02)

  tini, tend = 0.0, 4.

  tstep = 1/16

  tt = tini:tstep:tend

  ξ = set_variables("ξₓ ξᵥ", order=2, numvars=length(q0))

  tTM1 = tv["$(δq01)-$(δq02)"]

  xTM1 = tmv["$(δq01)-$(δq02)"]

  q0ξ = interval_rand(δq0)

  print("Initial conditions x0 = $(q0ξ + q0)")

  p = plot()

  for i in 2:length(tTM1)

    t0 = tTM1[i-1]

    trange = range(0, xTM1[1, i].dom.hi, length=20)

    fT = xTM1[1, i]

    lo, hi = fT.rem.lo, fT.rem.hi

    mid_y = [mid(fT(t)(symIbox)) for t in trange]

    upper_y = ([fT(t)(symIbox).hi + hi for t in trange] - mid_y)

    lower_y = (mid_y - [fT(t)(symIbox).lo + abs(lo) for t in trange])

    plot!(trange.+t0, mid_y, ribbon=(lower_y, upper_y), 

      c=:dimgray, label="", lw=3)

  end

  plot!(tt, t -> x(t, q0ξ .+ q0, g), label="Exact solution", lw=2, c=:red2)

  xlabel!("t")

  ylabel!("x(t)")

  #plot!(tt, t -> x(t, q0 .+ mid.(δq0) .+ sup.(δq0), g), color=:green,

  #  label="Sol inf")

  #plot!(tt, t -> x(t, q0 .+ mid.(δq0) .+ inf.(δq0), g), color=:red,

  #  label="Sol sup")

  p

end

References

[1] Florian Bünger, A Taylor model toolbox for solving ODEs implemented in MATLAB/INTLAB, Journal of Computational and Applied Mathematics (2020)

[2] Martin Berz & Kyoko Makino, Verified Integration of ODEs and Flows Using Differential Algebraic Methods on High-Order Taylor Models, Reliable Computing (1998)