# (Dis)Advantages of quintic (fifth order) vs. cubic (third order) Hermite interpolation?

**URL:** <https://discourse.julialang.org/t/dis-advantages-of-quintic-fifth-order-vs-cubic-third-order-hermite-interpolation/120383>\
**Category:** Numerics\
**Tags:** ode, interpolations\
**Created:** [September 28, 2024, 3:03pm UTC](https://discourse.julialang.org/t/dis-advantages-of-quintic-fifth-order-vs-cubic-third-order-hermite-interpolation/120383 "2024-09-28T15:03:54Z")\
**Posts on this page:** 19\
**Page:** 1

<div class="post-metadata">

**Author:** ![Marco\_Masi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marco_masi/32/212305_2.png) [@Marco\_Masi](https://discourse.julialang.org/u/Marco_Masi)\
**Post date:** [September 28, 2024, 3:03pm UTC](https://discourse.julialang.org/t/dis-advantages-of-quintic-fifth-order-vs-cubic-third-order-hermite-interpolation/120383/1 "2024-09-28T15:03:54Z")

</div>

After solving an ordinary differential equation (ODE) of a dynamical system and obtaining the solution in time with its first and second derivatives, I fit the data points with a QuinticHermiteSpline(ddu, du, u, t). It appears to be the most appropriate interpolation because it uses all the available information matching all the points and its derivatives. It is a spline where each piece is a fifth-degree polynomial. However, as far as I understand it, higher order interpolations aren’t always the best choice, because they could have large oscillations. Since I need to find for the most precise solution and make the most precise interpolations, I wonder if and how going for a fifth-degree polynomial is safe? Comparing the two interpolations, both seem to work well fitting the expected solution of the ODE, but they are (obviously) different on smaller scales. I’m not sure whether using a fifth-order is a more accurate fit vs. a third-order cubic Hermite interpolation, which, however, uses only the first derivative. Any idea when and where which of the two cases is preferable?

---

<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:** [September 28, 2024, 5:16pm UTC](https://discourse.julialang.org/t/dis-advantages-of-quintic-fifth-order-vs-cubic-third-order-hermite-interpolation/120383/2 "2024-09-28T17:16:30Z")

</div>

It’s fine. It can have some stability issues, but the main reason the ODE solvers use a cubic interpolation by default as a fallback is because it’s free, so 0 extra computations are required to construct it. A 5th order polynomial requires a bit more work.

But of course, neither of those are a good idea because you can construct higher order polynomials which are more stable with less work using the tableau values, which is why 5th order quintic is not used anywhere.

---

<div class="post-metadata">

**Author:** ![Marco\_Masi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marco_masi/32/212305_2.png) [@Marco\_Masi](https://discourse.julialang.org/u/Marco_Masi)\
**Post date:** [September 28, 2024, 8:07pm UTC](https://discourse.julialang.org/t/dis-advantages-of-quintic-fifth-order-vs-cubic-third-order-hermite-interpolation/120383/3 "2024-09-28T20:07:29Z")

</div>

Actually, I found this method so interesting because it requires no work at all. The methods illustrated [here](https://docs.sciml.ai/DataInterpolations/stable/methods/) are quite straightforward. Which higher order polynomials that are more stable are you thinking of in particular? Do they exactly interpolate the first and second derivative as well?

---

<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:** [September 28, 2024, 8:31pm UTC](https://discourse.julialang.org/t/dis-advantages-of-quintic-fifth-order-vs-cubic-third-order-hermite-interpolation/120383/4 "2024-09-28T20:31:26Z")

</div>

@Marco_Masi the piece you’re missing is that an ODE solver has more information than just the points. You also know the derivatives since the ODE formulation is `du/dt = f(u,p,t)` so if you have computed `f(u,p,t_0)` you know your derivative at `t_0`.

---

<div class="post-metadata">

**Author:** ![Marco\_Masi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marco_masi/32/212305_2.png) [@Marco\_Masi](https://discourse.julialang.org/u/Marco_Masi)\
**Post date:** [September 28, 2024, 9:00pm UTC](https://discourse.julialang.org/t/dis-advantages-of-quintic-fifth-order-vs-cubic-third-order-hermite-interpolation/120383/5 "2024-09-28T21:00:10Z")

</div>

I know. That’s why I’m asking. I’m using it for the interpolation.

---

<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:** [September 28, 2024, 9:53pm UTC](https://discourse.julialang.org/t/dis-advantages-of-quintic-fifth-order-vs-cubic-third-order-hermite-interpolation/120383/6 "2024-09-28T21:53:37Z")

</div>

> [@Marco\_Masi](#):
>
> The methods illustrated [here](https://docs.sciml.ai/DataInterpolations/stable/methods/) are quite straightforward.

They are straightforward but they require more computational work, i.e. more `f` evaluations.

> [@Marco\_Masi](#):
>
> Which higher order polynomials that are more stable are you thinking of in particular? Do they exactly interpolate the first and second derivative as well?

Yes they do, they are just ODE-solver specific. Note the ODE solver already has them setup, if you do `sol(t)` you’ll get an interpolation. In the print out of the solution it will tell you what accuracy it has.

---

<div class="post-metadata">

**Author:** ![Marco\_Masi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marco_masi/32/212305_2.png) [@Marco\_Masi](https://discourse.julialang.org/u/Marco_Masi)\
**Post date:** [September 29, 2024, 9:30am UTC](https://discourse.julialang.org/t/dis-advantages-of-quintic-fifth-order-vs-cubic-third-order-hermite-interpolation/120383/7 "2024-09-29T09:30:07Z")

</div>

Yes, but I made some tests and the ODE interpolation has the accuracy of the order of the cubicHermite (it looks like to be the very same.) While, in the case one integrates a second order ODE and makes the little effort to calculate the second derivative and interpolates with the quinticHermite, one gets a much more precise fit, about 4-5 orders of magnitude better!

---

<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:** [September 29, 2024, 9:56am UTC](https://discourse.julialang.org/t/dis-advantages-of-quintic-fifth-order-vs-cubic-third-order-hermite-interpolation/120383/8 "2024-09-29T09:56:59Z")

</div>

It depends on the method. Which second order ODE integrator? If you’re using a Verner method like Vern7 then it has some special interpolation, DPRKN as well.

---

<div class="post-metadata">

**Author:** ![Marco\_Masi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marco_masi/32/212305_2.png) [@Marco\_Masi](https://discourse.julialang.org/u/Marco_Masi)\
**Post date:** [September 29, 2024, 2:45pm UTC](https://discourse.julialang.org/t/dis-advantages-of-quintic-fifth-order-vs-cubic-third-order-hermite-interpolation/120383/9 "2024-09-29T14:45:34Z")

</div>

I tested Vern7 and vern9 vs. the analytic solution (comparison made with same tolerance and precision). Indeed, Vern7 is accurate as the QuinticHermite. Interestingly, however, Vern9 is less precise than Vern7. At least with my problem (non-linear ODE of the driven damped pendulum.) Any idea why? I would expect a higher order solver be more accurate… However, one can recover the same accuracy of Ver7 by resorting to QuinticHermite. I also tested Feagin14. It is less accurate than Vern7/9 but, again, one can recover the same accuracy by resorting to the QuinticHermite.  
The graphs below show the logarithmic divergence between the “real” solution (the exact analytic solution of the non-driven pendulum, it is an exponentially decaying function) and the numerical one. (Note how the blue lines are hidden behind the red ones, which means that the CubicHremite matches exactly the interpolation of the solvers. I suspect that internally that’s what they are using.)  
So, one might be tempted to go for Vern7. However, it is very slow compared to Faegin14 (which is almost ten times faster.) Thus, I might use the latter and then interpolate with Hermite.  
Other solvers that could do even better? Looking up the DifferentialEquations.jl there are far too many to test them all. Any thoughts?

![Test_Vern7_t10_tol30_prec3322](https://global.discourse-cdn.com/julialang/original/3X/1/a/1a8d89f825f6390369fcdd476da233e42ca57523.png)

![Test_Vern9_t10_tol30_prec3322](https://global.discourse-cdn.com/julialang/original/3X/4/e/4e958573b398223d27dde22f85d0b185dc5f7e46.png)

![Test_Feagin14_t10_tol30_prec3322](https://global.discourse-cdn.com/julialang/original/3X/7/f/7f26fc053b1997b37a7cb8e21d558af0bc70e3c0.png)

---

<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:** [September 29, 2024, 4:37pm UTC](https://discourse.julialang.org/t/dis-advantages-of-quintic-fifth-order-vs-cubic-third-order-hermite-interpolation/120383/10 "2024-09-29T16:37:56Z")

</div>

16 digits of accuracy is the floating point limit. You’re basically just hitting floating point truncation error in that plot

---

<div class="post-metadata">

**Author:** ![Marco\_Masi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marco_masi/32/212305_2.png) [@Marco\_Masi](https://discourse.julialang.org/u/Marco_Masi)\
**Post date:** [September 29, 2024, 6:56pm UTC](https://discourse.julialang.org/t/dis-advantages-of-quintic-fifth-order-vs-cubic-third-order-hermite-interpolation/120383/11 "2024-09-29T18:56:43Z")

</div>

Good point. But how could that be? I’m working with BigFloats and setprecision(3300) which, as I understand it, should be more than 1000 digits.

---

<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:** [September 29, 2024, 7:34pm UTC](https://discourse.julialang.org/t/dis-advantages-of-quintic-fifth-order-vs-cubic-third-order-hermite-interpolation/120383/12 "2024-09-29T19:34:02Z")

</div>

many of the ode solvers only have coefficients that satisfy the order conditions up to floating point tolerance.

---

<div class="post-metadata">

**Author:** ![Marco\_Masi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marco_masi/32/212305_2.png) [@Marco\_Masi](https://discourse.julialang.org/u/Marco_Masi)\
**Post date:** [September 29, 2024, 7:39pm UTC](https://discourse.julialang.org/t/dis-advantages-of-quintic-fifth-order-vs-cubic-third-order-hermite-interpolation/120383/13 "2024-09-29T19:39:43Z")

</div>

This is the difference between two solutions made with the same solver with the same tolerance but with different precisions. How can it appreciate differences smaller than 1e-16?

 ![forum](https://global.discourse-cdn.com/julialang/original/3X/e/4/e4fd8d2499ecda04158cc70110515a5d94e6550c.jpeg)

---

<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:** [September 29, 2024, 11:22pm UTC](https://discourse.julialang.org/t/dis-advantages-of-quintic-fifth-order-vs-cubic-third-order-hermite-interpolation/120383/14 "2024-09-29T23:22:17Z")

</div>

If you’re using bigfloats then yeah it can get that accurate. Though your plot is a bit odd: how could the interpolation with quintic give better results than the step points itself? The accuracy of the points shouldn’t change.

---

<div class="post-metadata">

**Author:** ![Marco\_Masi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marco_masi/32/212305_2.png) [@Marco\_Masi](https://discourse.julialang.org/u/Marco_Masi)\
**Post date:** [September 30, 2024, 8:22am UTC](https://discourse.julialang.org/t/dis-advantages-of-quintic-fifth-order-vs-cubic-third-order-hermite-interpolation/120383/15 "2024-09-30T08:22:43Z")

</div>

Oopss … my previous graph was wrong. It used different points in time than sol.t. Using the same time-vector, all converge on the same (log)error.

![Julia_discourse](https://global.discourse-cdn.com/julialang/original/3X/5/c/5cb68eb1fbd5a8051e8f509141528ccb48609330.jpeg)

However, the double precision limitation remains. It seems that passing BigFloat variables with no matter how many digits precision to the solvers and interpolators doesn’t make them work internally with the same typeof. That’s very unfortunate! Because I desperately need highest multi-precision solutions of non-linear ODEs which are highly sensitive to initial conditions (ergo, to “noisy” numerical inaccuracy). Maybe I’m doing something wrong (in the code or due to some fundamentally flawed reasoning.) Here I add the code, in case someone likes to check this out…

```julia
using OrdinaryDiffEq, Plots, DataInterpolations
setprecision(3322); # 3322 --> ca. 1000 digits

InPos=BigFloat("-1"); InMom=BigFloat("5"); # Initial conditions
Tin=BigFloat("0"); Tfin=BigFloat("10"); # Time span

# ODE
function pend(dy, y, p, t)
    (;gamma, w0) = p
    dy[1] = y[2]
    dy[2] = -gamma*y[2] - w0^2*y[1] # y[2]=y' ; y[1]=y
end

prob = ODEProblem(pend, maxiters = 1e7, [InPos,InMom], (Tin,Tfin), (gamma=BigFloat("1"), w0=sqrt(BigFloat("37"))));
@time sol = solve(prob, Feagin14(), abstol=1e-30, reltol=1e-30);

t=sol.t; L=length(t);
Tin=t[1]; Tfin=t[end];
u=[sol.u[i][1] for i in 1:L];
du=[sol.u[i][2] for i in 1:L];
gamma=BigFloat("1"); w0=sqrt(BigFloat("37"));
ddu= -gamma*du - w0^2*u; # This obtains the second derivative via dy[2] above

# This is the analytic solution
f(t)=-(1/7)*exp(-t/2)*(7*cos(7*sqrt(3)*t/2) - 3*sqrt(3)*sin(7*sqrt(3)*t/2));
y=map(t->f(t),t);

H3=CubicHermiteSpline(du, u, t);
H5=QuinticHermiteSpline(ddu, du, u, t);

ylog1=broadcast(log10, abs.(sol.(t,idxs=1)-f.(t)));
ylog2=broadcast(log10, abs.(H3.(t)-f.(t)));
ylog3=broadcast(log10, abs.(H5.(t)-f.(t)));

plot(t,ylog1, label=["log10|sol(t)-f(t)|"])
plot!(t,ylog2, label=["log10|H3(t)-f(t)|"])
plot!(t,ylog3, label=["log10|H5(t)-f(t)|"])
xlabel!("Time")
ylabel!("Log10|sol(t)-f(t)|")

```

---

<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:** [September 30, 2024, 12:44pm UTC](https://discourse.julialang.org/t/dis-advantages-of-quintic-fifth-order-vs-cubic-third-order-hermite-interpolation/120383/16 "2024-09-30T12:44:49Z")

</div>

you’re only computing your analytical solution to float64 precision

---

<div class="post-metadata">

**Author:** ![Marco\_Masi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marco_masi/32/212305_2.png) [@Marco\_Masi](https://discourse.julialang.org/u/Marco_Masi)\
**Post date:** [September 30, 2024, 2:44pm UTC](https://discourse.julialang.org/t/dis-advantages-of-quintic-fifth-order-vs-cubic-third-order-hermite-interpolation/120383/17 "2024-09-30T14:44:07Z")

</div>

Oh… right. I fixed that, and now it looks much better (same result but with a -33 error.) This begins to make more sense. Thank you so far.

---

<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:** [September 30, 2024, 4:46pm UTC](https://discourse.julialang.org/t/dis-advantages-of-quintic-fifth-order-vs-cubic-third-order-hermite-interpolation/120383/18 "2024-09-30T16:46:27Z")

</div>

glad to help! One last thing to consider is that you probably should set your number tolerance to ~1.5x the bits of your solver tolerance. More numeric precision will just be extra slow for no reason.

---

<div class="post-metadata">

**Author:** ![Marco\_Masi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marco_masi/32/212305_2.png) [@Marco\_Masi](https://discourse.julialang.org/u/Marco_Masi)\
**Post date:** [September 30, 2024, 5:23pm UTC](https://discourse.julialang.org/t/dis-advantages-of-quintic-fifth-order-vs-cubic-third-order-hermite-interpolation/120383/19 "2024-09-30T17:23:51Z")

</div>

Yes, I used 3300 bits only to be sure that the issue wasn’t related to a lack of precision. Making some further tests it appears that even less than 1/10 of that, would meet the requirement for 1e-30 tolerance. But I will continue to push the accuracy to the limits of my PC because ultimately I would like to check the scaling propagation of quantum noise in macroscopic systems, and being absolutely certain that it isn’t drowned in the numerical error propagation. Exploring the parameter space…