# DifferentialEquations extracting coefficients of dense interpolant

**URL:** <https://discourse.julialang.org/t/differentialequations-extracting-coefficients-of-dense-interpolant/93365>\
**Category:** Modelling & Simulations\
**Tags:** interpolations, differentialequation\
**Created:** [January 22, 2023, 3:54pm UTC](https://discourse.julialang.org/t/differentialequations-extracting-coefficients-of-dense-interpolant/93365 "2023-01-22T15:54:39Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![Bernard\_GODARD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bernard_godard/32/4155_2.png) [@Bernard\_GODARD](https://discourse.julialang.org/u/Bernard_GODARD)\
**Post date:** [January 22, 2023, 3:54pm UTC](https://discourse.julialang.org/t/differentialequations-extracting-coefficients-of-dense-interpolant/93365/1 "2023-01-22T15:54:39Z")

</div>

From the solution of a differential equation using Vern6, Vern7, Vern8 or Vern9, how can I extract the dense interpolant polynomial coefficients?

Of course I could evaluate the polynomial on a grid of size order+1 inside a step and then compute the interpolation polynomial using Lagrange’s method, but I do not want to evaluate the polynomial because the evaluation of that polynomial in Float64 arithmetics creates too much numerical noise for my application. I can solve this issue by either:

1. using BigFloat numeric type (tested in Julia)
2. evaluating the interpolation polynomial using compensated floating point arithmetics: [https://www-pequan.lip6.fr/~jmc/polycopies/Compensation-horner.pdf](https://www-pequan.lip6.fr/~jmc/polycopies/Compensation-horner.pdf) (I tested this successfully in a C++ code with Verner’s RK877 dense interpolant computed in double precision)

I would like to try method 2 in Julia but for that I need access to the interpolant coefficients.

Thank you.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [January 22, 2023, 4:23pm UTC](https://discourse.julialang.org/t/differentialequations-extracting-coefficients-of-dense-interpolant/93365/2 "2023-01-22T16:23:18Z")

</div>

> [@Bernard\_GODARD](#):
>
> From the solution of a differential equation using Vern6, Vern7, Vern8 or Vern9, how can I extract the dense interpolant polynomial coefficients?

`sol.k`.

> [@Bernard\_GODARD](#):
>
> 1. evaluating the interpolation polynomial using compensated floating point arithmetics: [https://www-pequan.lip6.fr/~jmc/polycopies/Compensation-horner.pdf](https://www-pequan.lip6.fr/~jmc/polycopies/Compensation-horner.pdf) (I tested this successfully in a C++ code with Verner’s RK877 dense interpolant computed in double precision)
> 
> I would like to try method 2 in Julia but for that I need access to the interpolant coefficients.

It would be easiest to just modify the code.

> <https://github.com/SciML/OrdinaryDiffEq.jl/blob/master/src/dense/interpolants.jl#L1659-L2603>

If it helps then I would definitely accept a PR to at least make it an option.

---

<div class="post-metadata">

**Author:** ![Bernard\_GODARD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bernard_godard/32/4155_2.png) [@Bernard\_GODARD](https://discourse.julialang.org/u/Bernard_GODARD)\
**Post date:** [January 30, 2023, 2:44pm UTC](https://discourse.julialang.org/t/differentialequations-extracting-coefficients-of-dense-interpolant/93365/3 "2023-01-30T14:44:46Z")

</div>

Thank you @ChrisRackauckas .

It looks to me that I just need to change Base.@evalpoly since this macro is only used in the interpolant (I want to leave the propagator steps unchanged. I understand the interpolants are used for event finding but I do not have events in this case).

```julia
$ find .julia/packages/OrdinaryDiffEq/ -type f -exec grep -q "evalpoly" {} \; -print
.julia/packages/OrdinaryDiffEq/W2xe3/src/dense/interpolants.jl

```

However it seems even if I redefine Base.@evalpoly before loading DifferentialEquations this has no effect on the @evalpoly used in DifferentialEquations.jl.

```julia
function TwoProductFMA(a::Float64,b::Float64)
    x = a * b;
    y = fma( a, b, -x )
    return x,y
end

function TwoSum( a::Float64, b::Float64 )
    x = a + b
    z = x - a
    y = ( a - ( x - z ) ) + ( b - z )
    return x,y
end

function chorner( x::Float64, coeffs)
    n = length(coeffs)
    r, corr = coeffs[n], zero(x)
    for i in 1:n-1
        r, a = TwoProductFMA( r, x );
        r, b = TwoSum( r, coeffs[n - i] );
        corr = corr * x;
        corr += a + b;
    end
    return r+corr
end

import Base.@evalpoly
macro evalpoly(z, p...)
    zesc, pesc = esc(z), esc.(p)
    :(chorner($zesc, ($(pesc...),)))
end

```

Is there a way to redefine Base.@evalpoly globally? Is there some compile cache that has to be manually cleanep up?

Note: I found in the Math.jl source there is already a compensated polynom evaluation called Base.exthorner but I cannot use it and I do not understand why:

```julia
julia> Base.exthorner
ERROR: UndefVarError: exthorner not defined

```

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [January 30, 2023, 3:06pm UTC](https://discourse.julialang.org/t/differentialequations-extracting-coefficients-of-dense-interpolant/93365/4 "2023-01-30T15:06:49Z")

</div>

It’s `Base.Math.exthorner`. Also it’s worth noting that it is only partially compensated in that it will get full accuracy when you have terms of decreasing magnitude (and possibly alternating signs) but will lose some precision if your polynomial is divergent. (This is because adding the other half of the compensation has a performance cost and `exthorner` is one I wrote for special function evaluation where the polynomials you are evaluating come from convergent series).

---

<div class="post-metadata">

**Author:** ![Bernard\_GODARD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bernard_godard/32/4155_2.png) [@Bernard\_GODARD](https://discourse.julialang.org/u/Bernard_GODARD)\
**Post date:** [January 30, 2023, 5:59pm UTC](https://discourse.julialang.org/t/differentialequations-extracting-coefficients-of-dense-interpolant/93365/5 "2023-01-30T17:59:29Z")

</div>

Redefining the macro in Base does not work but specializing Base.evalpoly works.

```julia

function TwoProductFMA(a::Float64,b::Float64)
    x = a * b;
    y = fma( a, b, -x )
    return x,y
end

function TwoSum( a::Float64, b::Float64 )
    x = a + b
    z = x - a
    y = ( a - ( x - z ) ) + ( b - z )
    return x,y
end;

function chorner( x::Float64, coeffs)
    n = length(coeffs)
    r, corr = Float64(coeffs[n]), zero(x)
    for i in 1:n-1
        r, a = TwoProductFMA( r, x );
        r, b = TwoSum( r, Float64(coeffs[n - i]) );
        corr = corr * x;
        corr += a + b;
    end
    return r+corr
end

import Base.evalpoly
evalpoly(x::Float64, p::Tuple)=chorner(x,p)
#evalpoly(x::Float64, p::AbstractVector) = chorner(x,p)

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [January 31, 2023, 7:52am UTC](https://discourse.julialang.org/t/differentialequations-extracting-coefficients-of-dense-interpolant/93365/6 "2023-01-31T07:52:49Z")

</div>

It seems the right thing to do might be to have a separate macro and an option to flip the choice?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 31, 2023, 1:38pm UTC](https://discourse.julialang.org/t/differentialequations-extracting-coefficients-of-dense-interpolant/93365/7 "2023-01-31T13:38:43Z")

</div>

> [@Bernard\_GODARD](#):
>
> I could evaluate the polynomial on a grid of size order+1 inside a step and then compute the interpolation polynomial using Lagrange’s method

If you are worried about numerical accuracy, you usually don’t want to use Lagrange’s formula, which is numerically unstable. Use e.g. barycentric interpolation, and where possible you want to avoid computing coefficients of polynomials in a monomial \{1,x,x^2,\ldots\} basis.

What are you trying to do?

---

<div class="post-metadata">

**Author:** ![Bernard\_GODARD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bernard_godard/32/4155_2.png) [@Bernard\_GODARD](https://discourse.julialang.org/u/Bernard_GODARD)\
**Post date:** [February 26, 2023, 5:15pm UTC](https://discourse.julialang.org/t/differentialequations-extracting-coefficients-of-dense-interpolant/93365/8 "2023-02-26T17:15:00Z")

</div>

> It seems the right thing to do might be to have a separate macro and an option to flip the choice?

That would certainly be nice. For my use case redefining Base.evalPoly worked.

---

<div class="post-metadata">

**Author:** ![Bernard\_GODARD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bernard_godard/32/4155_2.png) [@Bernard\_GODARD](https://discourse.julialang.org/u/Bernard_GODARD)\
**Post date:** [February 26, 2023, 6:41pm UTC](https://discourse.julialang.org/t/differentialequations-extracting-coefficients-of-dense-interpolant/93365/9 "2023-02-26T18:41:24Z")

</div>

@stevengj I do not want to use Lagrange Interpolation. I was just saying that I could obtain the polynomial coefficients using a Lagrange interpolator if I could evaluate the ODE solution interpolator to any accuracy (eg with BigFloat). But ideally I just want to get the interpolator coefficients. My goal was to assess different formulas to evaluate the interpolation polynomial (compensated arithmetic) to get rid of numerical noise in my model. The noise in my model is amplified by cancellation because the main term in the model is almost of the form of a finite difference:

\frac{f(u(t+T))-f(u(t-T))}{2T}

where u(t) is the ODE solution.

Ideally we get rid of the numerical noise by using a Taylor serie in T  
and that works well but requires accurate derivatives (automatic or analytic) for f \circ u and that is required in some case.

But already reducing the noise in the ODE solution interpolator polynomial evaluation helps a lot.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [February 26, 2023, 7:51pm UTC](https://discourse.julialang.org/t/differentialequations-extracting-coefficients-of-dense-interpolant/93365/10 "2023-02-26T19:51:16Z")

</div>

> [@Bernard\_GODARD](#):
>
> works well but requires accurate derivatives (automatic or analytic) for f \circ uf∘uf \circ u and that is required in some case.

The interpolations have the derivative forms written out as well with sol(t,Val{1})