Extrapolation Methods

Compute p(x) where p is the unique polynomial of degree length(xi), such that p(x[i]) = y[i]) for all i. Call with

aitken_neville!(x::AbstractVector, t::Real, ti::AbstractVector, xi::AbstractMatrix)

where

• x: evaluation value
• t: evaluation point
• ti: interpolation nodes
• xi: interpolation values

Euler extrapolation method with arbitrary order p.

Solves the ordinary differential equation

\[\begin{aligned} \dot{x} &= v(t, x) , & x(t_0) &= x_0 , \end{aligned}\]

for $x_1 = x(t_1)$, and is called with

extrapolate!(t₀, x₀, t₁, x₁, problem, EulerExtrapolation(s))

where

• t₀: initial time
• t₁: final time
• x₀: initial value $x_0 = x(t_0)$
• x₁: final value $x_1 = x(t_1)$
• problem: ODEProblem whose solution to extrapolate
• s: number of interpolations (order $p=s+1$)

Hermite's Interpolating Polynomials

Implements a two point Hermite inter-/extrapolation function which passes through the function and its first derivative for the interval $[0,1]$. The polynomial is determined by four constraint equations, matching the function and its derivative at the points $0$ and $1$.

Call with one of the following methods

extrapolate!(t₀, x₀, ẋ₀, t₁, x₁, ẋ₁, t, x, HermiteExtrapolation())
extrapolate!(t₀, x₀, ẋ₀, t₁, x₁, ẋ₁, t, x, ẋ, HermiteExtrapolation())
extrapolate!(t₀, x₀, t₁, x₁, t, x, v, HermiteExtrapolation())
extrapolate!(t₀, x₀, t₁, x₁, t, x, ẋ, v, HermiteExtrapolation())
extrapolate!(t₀, x₀, t₁, x₁, t, x, problem, HermiteExtrapolation())
extrapolate!(t₀, x₀, t₁, x₁, t, x, ẋ, problem, HermiteExtrapolation())

where

• t₀: first sample time $t_0$
• x₀: first solution value $x_0 = x(t_0)$
• ẋ₀: first vector field value $ẋ_0 = v(t_0, x(t_0))$
• t₁: second sample time $t_1$
• x₁: second solution value $x_1 = x(t_1)$
• ẋ₁: second vector field value $ẋ_1 = v(t_1, x(t_1))$
• t: time $t$ to extrapolate
• x: extrapolated solution value $x(t)$
• ẋ: extrapolated vector field value $ẋ(t)$
• v: function to compute vector field with signature v(ẋ,t,x)
• problem: ODEProblem whose vector field to use

Derivation

The interpolation works as follows: Start by defining the 3rd degree polynomial and its derivative by

\[\begin{aligned} g(x) &= a_0 + a_1 x + a_2 x^2 + a_3 x^3 , \\ g'(x) &= a_1 + 2 a_2 x + 3 a_3 x^2 , \end{aligned}\]

and apply the constraints

\[\begin{aligned} g(0) &= f_0 & & \Rightarrow & a_0 &= f_0 , \\ g(1) &= f_1 & & \Rightarrow & a_0 + a_1 + a_2 + a_3 &= f_1 , \\ g'(0) &= f'_0 & & \Rightarrow & a_1 &= f'_0 , \\ g'(1) &= f'_1 & & \Rightarrow & a_1 + 2 a_2 + 3 a_3 &= f'_1 . \\ \end{aligned}\]

Solving for $a_0, a_1, a_2, a_3$ leads to

\[\begin{aligned} a_0 &= f_0 , & a_1 &= f'_0 , & a_2 &= - 3 f_0 + 3 f_1 - 2 f'_0 - f'_1 , & a_3 &= 2 f_0 - 2 f_1 + f'_0 + f'_1 , \end{aligned}\]

so that the polynomial $g(x)$ reads

\[g(x) = f_0 + f'_0 x + (- 3 f_0 + 3 f_1 - 2 f'_0 - f'_1) x^2 + (2 f_0 - 2 f_1 + f'_0 + f'_1) x^3 .\]

The function and derivative values can be factored out, so that $g(x)$ can be rewritten as

\[g(x) = f_0 (1 - 3 x^2 + 2 x^3) + f_1 (3 x^2 - 2 x^3) + f'_0 (x - 2 x^2 + x^3) + f'_1 (- x^2 + x^3) ,\]

or in generic form as

\[g(x) = f_0 a_0(x) + f_1 a_1(x) + f'_0 b_0(x) + f'_1 b_1(x) ,\]

with basis functions

\[\begin{aligned} a_0 (x) &= 1 - 3 x^2 + 2 x^3 , & b_0 (x) &= x - 2 x^2 + x^3 , \\ a_1 (x) &= 3 x^2 - 2 x^3 , & b_1 (x) &= - x^2 + x^3 . \end{aligned}\]

The derivative $g'(x)$ accordingly reads

\[g'(x) = f_0 a'_0(x) + f_1 a'_1(x) + f'_0 b'_0(x) + f'_1 b'_1(x) ,\]

with

\[\begin{aligned} a'_0 (x) &= - 6 x + 6 x^2 , & b'_0 (x) &= 1 - 4 x + 3 x^2 , \\ a'_1 (x) &= 6 x - 6 x^2 , & b'_1 (x) &= - 2 x + 3 x^2 . \end{aligned}\]

The basis functions $a_0$and $a_1$ are associated with the function values at $x_0$ and $x_1$, respectively, while the basis functions $b_0$ and $b_1$ are associated with the derivative values at $x_0$ and $x_1$. The basis functions satisfy the following relations,

\[\begin{aligned} a_i (x_j) &= \delta_{ij} , & b_i (x_j) &= 0 , & a'_i (x_j) &= 0 , & b'_i (x_j) &= \delta_{ij} , & i,j &= 0, 1 , \end{aligned}\]

where $\delta_{ij}$ denotes the Kronecker-delta, so that

\[\begin{aligned} g(0) &= f_0 , & g(1) &= f_1 , & g'(0) &= f'_0 , & g'(1) &= f'_1 . \end{aligned}\]

Midpoint extrapolation method with arbitrary order p.

For an ODEProblem, this solves the ordinary differential equation

\[\begin{aligned} \dot{x} &= v(t, x) , & x(t_0) &= x_0 , \end{aligned}\]

for $x_1 = x(t_1)$, and is called with

extrapolate!(t₀, x₀, t₁, x₁, ::ODEProblem, MidpointExtrapolation(s))

where

• t₀: initial time
• x₀: initial value $x_0 = x(t_0)$
• t₁: final time
• x₁: final value $x_1 = x(t_1)$
• s: number of interpolations (order $p=2s+2$)

For a PODEProblem or HODEProblem, this solves the partitioned ordinary differential equation

\[\begin{aligned} \dot{q} &= v(t, q, p) , & q(t_0) &= q_0 , \\ \dot{p} &= f(t, q, p) , & p(t_0) &= p_0 , \end{aligned}\]

for $q_1 = q(t_1)$ and $p_1 = p(t_1)$m and is called with

extrapolate!(t₀, q₀, p₀, t₁, q₁, p₁, ::PODEProblem, MidpointExtrapolation(s))
extrapolate!(t₀, q₀, p₀, t₁, q₁, p₁, ::HODEProblem, MidpointExtrapolation(s))

where

• t₀: initial time
• q₀: initial position $q_0 = q(t_0)$
• p₀: initial momentum $p_0 = p(t_0)$
• t₁: final time
• q₁: final position $q_1 = q(t_1)$
• p₁: final momentum $p_1 = p(t_1)$
• s: number of interpolations (order $p=2s+2$)

Similarly, for a IODEProblem or LODEProblem, this solves the explicit dynamical equation

\[\begin{aligned} \dot{q} &= v(t, q) , & q(t_0) &= q_0 , \\ \dot{p} &= f(t, q, v) , & p(t_0) &= p_0 , \end{aligned}\]

corresponding to the implicit problem, for $q_1 = q(t_1)$ and $p_1 = p(t_1)$, and is called with

extrapolate!(t₀, q₀, p₀, t₁, q₁, p₁, ::IODEProblem, MidpointExtrapolation(s))
extrapolate!(t₀, q₀, p₀, t₁, q₁, p₁, ::LODEProblem, MidpointExtrapolation(s))

where

• t₀: initial time
• q₀: initial position $q_0 = q(t_0)$
• p₀: initial momentum $p_0 = p(t_0)$
• t₁: final time
• q₁: final position $q_1 = q(t_1)$
• p₁: final momentum $p_1 = p(t_1)$
• s: number of interpolations (order $p=2s+2$)