[Math] Determining the parameters of a differential equation

Question

Let's assume:

$$ay''(t)+by'(t)+cy(t)=f(t)$$

And we are given $n$ points $\{(t_1,y_1),(t_2,y_2),\ldots,(t_n,y_n)\}$, where $y_i=y(t_i)+\epsilon_i$.

By approximation of the differential equation, I mean finding some unknown parameters of the ordinary differential equation (ODE) itself and not the solution of the ODE.

My goal is to approximate $a,b,c$, such that $$\sum_i \|y(t_i)-y_i\|_2^2$$ is minimized.

Note 1: We might have extra constrains/approximations for $y''$ and $y'$. But those can be ignored in the first place.

Note 2: The objective $\sum_i \|y(t_i)-y_i\|_2^2$ is just a tentative one, if there is any other criterion that results in some sort approximation of the ODE that would be okay as well.

Note 3: It is okay to assume $\epsilon_i$'s are all equal to zero, if this assumption simplifies the problem.

I tried to get some ideas by looking into the literature of stochastic differential equations (SDEs), but it seems like in the SDE's literature the focus is on the noisiness of the parameters. (Which I'm not sure if it's the case here.)

I know we may be able to solve the differential equation parametrically (i.e. finding $y(t)$ explicitly) and then approximate the parameters by numerically minimizing the summation: $\sum_i \|y(t_i)-y_i\|_2^2$. However, I'm trying not to do that.

Thanks for the help!

Answer 1

Added 13 Oct, better than the below:

Derivative estimates are very sensitive to noise in the data, because high frequencies are amplified. Non-uniform spacing makes this even worse -- dividing by $t_{i+1} - t_i \approx 0$ is a bad thing to do. So don't do that then.
Instead, fit a spline to the data, a piecewise cubic with say 5 or 10 pieces. Then on each piece,
$\qquad spline(t) = a_0 + a_1 t + a_2 t^2 + a_3 t^3 $
$\qquad spline'(t) = \qquad a_1 + 2 a_2 t + 3 a_3 t^2 $
$\qquad spline''(t) = \qquad\qquad 2 a_2 + 6 a_3 t $

Sample these on whatever grid you like -- the origial $t_i$, or a uniform grid. (Bspline programs like scipy.interpol.splev will do $spline'$ and $spline''$ too.) Then you could fit $a\ b\ c$ with least squares. Or, Harmonic inversion looks nice (but I haven't used it):

the FDM algorithm works by considering the time-series to be the autocorrelation function of a fictitious dynamical system ... For $M$ data points and $J$ frequencies, the time required is $O(M J) + O(J^3)$

---

First do the simpler case of $t_i = 1\ 2\ 3 ...$ equally spaced at integers.
Notation: call $y(t_i)\ y_i$, and the given data values $Y_i$ .
Approximate the derivatives with Finite differences:
$\qquad y''_i \approx (y_{i+1} - 2 y_i + y_{i-1}) / 2 $
$\qquad y'_i \approx y_{i+1} - y_i \ \ $ or $(y_{i+1} - y_{i-1}) / 2 $
Then fitting the data
$\qquad Y_i \approx a\ (y_{i+1} - 2 y_i + y_{i-1}) / 2 \ + b\ (y_{i+1} - 2 y_i + y_{i-1}) / 2 \ + c\ y_i - f_i $
is a Least_squares problem: set up an $n \times 3$ matrix $A$, and fit
$\qquad A x \approx Y - f,\ x = [a\ b\ c] $ .
(The first row of A needs $y_0 \approx y_1$, and similarly the last row needs $y_{n+1} \approx y_n $ .)

Least squares is worth knowing about, because

• there are solid methods and programs
• it's easy to add constraints
• one can see outliers, way-off $Y_i$, and perhaps trim them
• it gives error bars / confidence intervals for $a\ b\ c$
  (if errors are normally distributed -- plot them).

---

Corrected 16 October (thanks @Alt):
splines, above, are much less sensitive to $t_{i+1} - t_i \approx 0$ than using just 3 or 5 neighbours. But if you must, here's how, in Python not equations. The gist has a little test case.

""" approximate derivatives for non-uniformly spaced data
    deriv1, deriv2( x, t ) with t0 < t1 < ...
    Very sensitive to noise, or delta-t near 0
    Spline, then take derivatives, is much better; see scipy.interpolate.splev( der= )
"""
# http://google.com/search?q="finite-difference non-uniform site:stackexchange.com"

from __future__ import division
import numpy as np

__version__ = "2015-10-18 oct  gist.github.com/denis-bz"

#...............................................................................
def deriv1( x, t ):
    return np.gradient(x) / np.gradient(t)  # (x+ - x-) / (t+ - t-)
    # return (ldiff(x) / ldiff(t)  +  rdiff(x) / rdiff(t)) / 2  ?

def deriv2( x, t ):
        # way off where delta-t is near 0
        # 5 neighbours, not 3
    return deriv1( deriv1( x, t ), t )