Approximate derivative within ODE function.

Hi, I am solving a pretty big system of differential equations using 'ode15s'. Within the ODE function, I would like to calculate the second derivative of a variable for which I don't have an analytical expression (even though I know its first derivative). I wonder if there's any possible way to calculate an approximate second derivative for that variable, considering that this last appears in other model's equations.

Thanks in advance, Davide

[Star Strider]

Davide Masiello’s follow-up Comment to my (now deleted) Answer:

I need to calculate the second derivative in the ODE function script because it appears as a member of other equations constituting the model, as in the following example:

where the second derivative of x is not given to know analytically. In the past, I used to stop the integration at a fixed and arbitrarily small time step to calculate the approximate derivative and then keep integrating with the last calculated value as new initial one, but I'd like to avoid this as the code runs correctly only with a varying step size. That's why I am wondering whether there's a more elegant solution to this problem.

Answer 1

[Jan]

If you have the first derivative, you can approximate the second derivative by the one or two sided quotient of differences:

ddf = (df(x+h) - df(x))   / h
ddf = (df(x+h) - df(x-h)) / (2*h)

The 2nd is more precise and usually twice as expensive assuming that df(x) is calculated at all.

A serious problem is the choice of a suitable small step h. It must be small enough to avoid that the 3rd derivative of the function changes the value of the 2nd derivative too much. But if it is too small, the cancellation problem occurs: The difference df(x+h)-df(x) has to few significant digits and dividing by a tiny h amplifies the rounding errors. As a rule of thumb df(x+h)-df(x) should be in the magnitude of sqrt(eps). An optimal h can be found in an iterative process only. For this using the former value of h is a good start point.

If the function to be integrated is a vector, you need different values for h for each component.

The problem of automatic differentiation is not trivial and you can write a PhD thesis about this topic. For details I suggest: Gill, Murray, Wright: Practical Optimization, 1981, Chapter 8.6.

[Davide Masiello]

Hi Jan, thanks for your exhaustive answer. It looks that would help me solving the problem, but in my particular case the step of integration "h" is decided by the Matlab internal routine of integration. Would you still suggest to find an optimal h value by trial and error or, at your knowledge, there's a way to get the integration step of the ode15s routine while it runs?

[Jan]

There are different "h" steps: The stepsize required to calculate the derivative accurately is not the one of the integration. In the given reference the method to determine the optimal stepsize for the determination of the derivative is explained exhaustively. It is not "trial and error".