I would like to compare my simulations in Mathematica to measured data. I measure acceleration but for simulations the starting point is displacement. How does one get accurate simulations of higher derivatives of the solutions to differential equations?

Here is a simple example of a simulation where we have the exact solution. (My actual problems are much more complicated). First I define a few parameters and then get exact results for displacement, velocity and acceleration.

```
vals = {
a -> 100, (* frequency *)
b -> 10, (* damping *)
y0 -> 1, (* initial displacement *)
v0 -> 0 (* initial velocity *)
};
tmax = 0.5 ;(* simulation time *)
SetOptions(Plot, PlotRange -> All, ImageSize -> 2.5 72);
d0 = y(t) /.
First@DSolve({y''(t) + 2 b y'(t) + (2 π a)^2 y(t) == 0,
y(0) == 1, y'(0) == 0}, y(t), t);
de = d0 /. vals;
ve = D(d0, t) /. vals;
ae = D(d0, {t, 2}) /. vals;
Row({
Plot(de, {t, 0, tmax}),
Plot(ve, {t, 0, tmax}),
Plot(ae, {t, 0, tmax})
})
```

The three plots are displacement, velocity and acceleration. The scales are very different but are noted in the next stage.

**One equation**

I use one equation in the simulation using `NDSolve`

and look at the errors between the displacement velocity and acceleration. The errors are normalised so that the functions being differenced are of order 1.

```
vmax = 600; amax = 400000;
d1 = y /. First@NDSolve({
y''(t) + 2 b y'(t) + (2 (Pi) a)^2 y(t) == 0,
y(0) == y0, y'(0) == v0} /. vals, y, {t, 0, tmax});
Row({Plot(d1(t) - de, {t, 0, tmax}),
Plot((d1'(t) - ve)/vmax, {t, 0, tmax}),
Plot((d1''(t) - ae)/amax, {t, 0, 0.5})})
```

The errors in displacement and velocity are of order 10^-7 which is reasonable but the error in acceleration is of order 10^-5. The errors in acceleration are particularly uneven with significant spikes.

**Two equations**

The differential equation is now split into two equations, one for displacement and one for velcity.

```
{d2, v2} = {y, v} /. First@NDSolve({
y'(t) == v(t),
v'(t) + 2 b v(t) + (2 (Pi) a)^2 y(t) == 0,
y(0) == y0, v(0) == v0} /. vals, {y, v}, {t, 0, tmax});
Row({
Plot(d2(t) - de, {t, 0, tmax}),
Plot((v2(t) - ve)/vmax, {t, 0, tmax}),
Plot((v2'(t) - ae)/amax, {t, 0, tmax})
})
```

The accuracy is worse for the displacement and velocity and no better for the acceleration.

**Three equations**

Here we have differential equations for the velocity and acceleration. The interaction between these derivatives is expressed as an algebraic equation.

```
{d3, v3, ac3} = {y, v, ac} /. First@NDSolve({
v'(t) == ac(t),
y'(t) == v(t),
ac(t) + 2 b v(t) + (2 π a)^2 y(t) == 0,
y(0) == y0, v(0) == v0} /. vals, {y, v, ac}, {t, 0, tmax});
Row({
Plot(d3(t) - de, {t, 0, tmax}),
Plot((v3(t) - ve)/vmax, {t, 0, tmax}),
Plot((ac3(t) - ae)/amax, {t, 0, tmax})
})
```

For this simulation the three results are of order 10^-6 and there are no spikes on any of the simulations. I guess the issues I am seeing are to due to the type of interpolation used. Derivatives of interpolation functions do not work well.

This suggests that the third approach is best if I want about equal errors for displacement, velocity and acceleration. Is this the best approach? Are there other approaches to getting a good second derivative?

If one had a third order, or forth order differential equation (I sometimes do) should one extend this method?