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?
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.
Actually, I found this method so interesting because it requires no work at all. The methods illustrated here 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?
@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.
I know. That’s why I’m asking. I’m using it for the interpolation.
They are straightforward but they require more computational work, i.e. more f evaluations.
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.
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!
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.
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?
16 digits of accuracy is the floating point limit. You’re basically just hitting floating point truncation error in that plot
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.
many of the ode solvers only have coefficients that satisfy the order conditions up to floating point tolerance.
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?
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.
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.
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…
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)|")
you’re only computing your analytical solution to float64 precision
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.
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.
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…