4
$\begingroup$

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}]
  }]

Mathematica graphics

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, tmax}]}]

Mathematica graphics

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}]
  }]

Mathematica graphics

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}]
  }]

Mathematica graphics

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?

Edit

Following the suggestion using the answer from I M I have investigated setting the interpolation order. The information from Help is not very clear it states

Use InterpolationOrder->All to get interpolation the same order as the method:

Can anyone explain what this means? What are orders of the method?

I have rerun the case with one equation. This gives:

vmax = 600; amax = 400000;
d1a = 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}, 
     InterpolationOrder -> All];

Row[{Plot[d1a[t] - de, {t, 0, tmax}],
  Plot[(d1a'[t] - ve)/vmax, {t, 0, tmax}],
  Plot[(d1a''[t] - ae)/amax, {t, 0, tmax}]}]

Mathematica graphics This is now good with similar errors for all derivatives.

Edit 2 Here is another approach suggested by Michael E2 that uses InterpolationOrder-> All, Method->"ImplicitRungeKutta"

vmax = 600; amax = 400000;
d1b = 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}, 
     InterpolationOrder -> All, Method -> "ImplicitRungeKutta"];

Row[{Plot[d1b[t] - de, {t, 0, tmax}],
  Plot[(d1b'[t] - ve)/vmax, {t, 0, tmax}],
  Plot[(d1b''[t] - ae)/amax, {t, 0, tmax}]}]

Mathematica graphics This gives an error of 10^-15 which is as much as one can hope for.

A further suggestion by J.M's ttechnical difficulties is to use Interpolation -> All, Method->"StiffnessSwitching"

vmax = 600; amax = 400000;
d1c = 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}, 
     InterpolationOrder -> All, Method -> "StiffnessSwitching"];

Row[{Plot[d1c[t] - de, {t, 0, tmax}],
  Plot[(d1c'[t] - ve)/vmax, {t, 0, tmax}],
  Plot[(d1c''[t] - ae)/amax, {t, 0, tmax}]}]

Mathematica graphics This is rather disappointing as there are spikes showing a huge error. However, if you zoom in, ignoring the spikes the error is about 10^-8.

$\endgroup$
1
$\begingroup$

You can increase interpolation order in NDSolve[], not sure, looks like it should be order + 2 or use All:

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},InterpolationOrder -> 4 (*All*)];
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}]}]

Edit

Thanks to @xzczd's comment I've remembered, "InterpolationOrder" is related to "DifferenceOrder", not ODE order.

You can use "ImplicitRungeKutta" or "ExplicitRungeKutta" methods, with custom "DifferenceOrder" and "InterpolationOrder" -> All.

For example, this is order 20 adaptive step implicit RK (also symplectic):

Method->{"ImplicitRungeKutta",DifferenceOrder->20,"Coefficients"->"ImplicitRungeKuttaGaussCoefficients","ImplicitSolver"->{"Newton",AccuracyGoal->MachinePrecision,PrecisionGoal->MachinePrecision,"IterationSafetyFactor"->1}}

You can also use fixed step method:

Method->{"TimeIntegration"->{"FixedStep",Method->{"ImplicitRungeKutta",DifferenceOrder->10,"Coefficients"->"ImplicitRungeKuttaGaussCoefficients","ImplicitSolver"->{"Newton",AccuracyGoal->MachinePrecision,PrecisionGoal->MachinePrecision,"IterationSafetyFactor"->1}}}}

For fixed step also set "StartingStepSize", "MaxStepFraction" and "MaxSteps" options.

| improve this answer | |
$\endgroup$
  • $\begingroup$ Thanks I have updated the question showing the effect of your suggestion. $\endgroup$ – Hugh Jul 8 at 14:40

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.