2
$\begingroup$

I am going to calculate a trajectory, by using a pre-calculated vector field. The values of the field are known on a grid which is quadratic in the horizontal direction, un-evenly spaced in the vertical, and evenly spaced in time.

I'm planning to use a 4th order Runge-Kutta integrator, and my question is if I should simply select my time step such that the three different times in the Runge-Kutta calculation exactly matches three consecutive time steps in the data, or if I will gain anything by implementing some kind of interpolation scheme to allow using a shorter time step in the integrator.

The spatial and temporal resolution of the pre-calculated dataset is such that the trajectory will seldom move further than the spatial step during one time step.

|cite|improve this question
$\endgroup$
2
$\begingroup$

Basically what RK does can be seen as some sort of interpolation. As you can see on Figure 1, the advantage of RK compared to Euler is its greater stability region so I think you should run a series of tests to check whether you are indeed inside the stability region. If you chose to interpolate first on your own, then you are going to introduce additional correlations that should not be there... I think it would be better to reduce the timestep in the computation of the vector field if you can.

|cite|improve this answer
$\endgroup$
  • $\begingroup$ I thought about testing, and I probably will, but it would be nice to be have a more general proof or argument. Recalculating the input data is not an option, I'm afraid, as that would take several months. $\endgroup$ – Tor Mar 20 '16 at 15:36
1
$\begingroup$

It would likely be best to not use the interpolation since I assume you'd default to cubic splines which is $\mathcal{O}(h^3)$ whereas RK4 is $\mathcal{O}(h^4)$. You'd have to do an error analysis to see if this would decrease your accuracy, and I think it's likely that it would. That said, if you do a really high order interpolation, then you're fine. To know how high of an order you need to do, you really have to look at the equation, put in the error term, check the expansion, and see where it ends up.

|cite|improve this answer
$\endgroup$
  • $\begingroup$ So, if I don't interpolate, that means I'm using the value at the nearest grid point, and the value changes discontinuously whenever I move to another gridpoint being the closest (also called nearest neighbour interpolation). Surely, that also introduces an error, and my gut feeling says this error is larger than that of an advanced interpolation scheme such as cubic splines. $\endgroup$ – Tor Mar 23 '16 at 20:01
  • $\begingroup$ I assumed that you had a way to ensure you'd step to the values in the grid, which would be errorless in function calculation but likely with a larger stepsize, meaning you'd have to look at the temporal error vs the grid error to make the determination. However, with your explanation in the comment, what you're actually saying is you're interpolating (via nearest neighbor interpolation). Yes this will induce error, and so your result will have an error term due to interpolation (I believe times $\mathcal{O}(h^4}$), check it) but the end result should be fine. $\endgroup$ – Chris Rackauckas Mar 23 '16 at 20:15
  • $\begingroup$ Though I would suggest using a better interpolation scheme like cubic splines. Nearest neighbor is discontinuous so your actual integration function is discontinuous, meaning all bets are off in terms of order of accuracy calculations for the temporal integration. Ideally you'd use a fourth order interpolation scheme here. $\endgroup$ – Chris Rackauckas Mar 23 '16 at 20:16
  • $\begingroup$ I see that I could probably have made the original question more clear. When I wrote it, I didn't think of nearest neighbour as interpolation. I thought about using quadrilinear interpolation, since then at least the values (but not their derivatives) will be continuous. Cubic splines would be nice, but whether it's viable depends entirely on performance. At each time step of my input data, I have two variables, each with about 250 x 150 x 45 grid points. The calculations I'm doing aren't that complicated, so I could easily see the splines calculation taking a significant fraction of the time. $\endgroup$ – Tor Mar 25 '16 at 8:09
  • $\begingroup$ At each timestep you can grab the cube that you're in and average between those points taking into account the distance for a linear interpolation, or double the cube for a local quadratic interpolation (you'd have to derive a formula here, I don't know where you'd find it, but it wouldn't be hard). $\endgroup$ – Chris Rackauckas Mar 25 '16 at 14:42
0
$\begingroup$

I did a bunch of tests, and it seems the answer to my question is that it does have an effect to use interpolation with Runge-Kutta. However, using a shorter timestep in the integrator than in the input data does not seem to have a large effect (I still need to interpolate, though, since the Runge-Kutta method evaluates the derivatives at the halfway point in time).

I have only tested this for a case where the ratio between spatial and temporal resolution in the input data are such that I will usually move only a fraction (1/10 - 1/5) of the spatial grid spacing during one timestep.

Interestingly, the choice of interpolation scheme seems to have an effect in some cases.

Particle transport 1

Particle transport 1

In the figures, I have transported particles for 2 days (top) and 3 days (bottom), using 4th order Runge-Kutta, with the interpolation scheme given by the legend. The x and y coordinates are in meters. For each scheme, I have used different timesteps: 2 hours, 1 hour, 0.5 hours and 0.25 hours. The timestep in the input data is 1 hour. Essentially, lowering the timestep below that of the input data has very little effect on the results. Going from linear interpolation to cubic splines has a large effect in one of the two cases shown below, though. It took about three times as long to run the program with splines interpolation.

I guess the only way to be sure which approach is "most correct" would be to get higher resolution data from the same source, and then compare trajectories derived from the high and low resolution data.

|cite|improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.