Solving Second Order Differential Equations with Discrete Time Data

Let's say I have a second order differential equation that is a function of discrete time series data that I have.

A simple example:

x** + x* = F(t) where, F(t) = rand(5000,1); <- some time series data sampled at a known frequency.

Does anyone know how I can solve this numerically for x ?

It seems that the examples of solving these equations all assume we know the actual function F(t), but in my case it is discrete data.

Thanks,

Kevin

Answer 1

[Walter Roberson]

Is that a system of equations,

 g''(t) + g'(t) = F(t)

i.e.,

 diff(g(x),x,x)(t) + diff(g(x),x)(t) = F(t)

If so then an equation is still needed to solve, and it seems to me that equation would need to have as high a order as the number of boundary conditions. Unless, that is, you want to be fitting the constants somehow ?

[Kevin]

Yes, that equation is a simple second order differential equation.

There must be a way to numerically solve this system by knowing the discrete time series data F(t) (but not its time dependant function).

[Walter Roberson]

If you have an unknown function of the form

diff(g(x), x, x)+diff(g(x), x) = C

with C an unknown constant, and you have values of the form

diff(g(x),x,x)(t) + diff(g(x),x)(t) = F(t), t = 1, 2, 3, ...

as a series of boundary conditions, then g(x) must be of the form

g(x) = -exp(-x)*C1 + C*x + C2

where C1 and C2 are constants of integration whose values are constrained by the boundary conditions.

You can see by examination that knowing any three of the boundary conditions would be enough to determine the unknown values C (the constant in the differential equation) and C1 and C2 (the constants of integration).

Thus, if you have more than 3 time points, your differential equation is over-determined.

For example, F(1)=2, F(2)=3, F(3)=5 is enough to uniquely identify

g(x) = (-exp(-x-1) + exp(-x+1) + exp(-x-2) - 1 + (3*exp(-3) - 5*exp(-2) + 2)*x) / (2*exp(-3) + 1 -3*exp(-2))

[Walter Roberson]

It is possible to find the polynomial f(x) of order N-1, N being the number of points in the time series, with f(1)=F(1), f(2)=F(2) and so on; this can be done through any of a number of techniques including constructing the coefficient matrix and using the backslash operator. This will, unfortunately, be amazingly inaccurate for higher N if you proceed numerically.

Once you have the polynomial and it is of the form [sum over K=1 to N of C(K)*x^(N-K)] where C denotes an array of constants, then it is possible to mechanically create the solution to

diff(g(x),x,x) + diff(g(x),x) = f(x)

as a sum of terms involving pochhamer() calculations and the constants. pochhammer(K,N) is gamma(K+N)/gamma(N) which is (K+N-1)!/(N-1)!

For example for N=8, one of the terms is

(1/5*(840*C(1)-210*C(2)+42*C(3)-7*C(4)+C(5)))*x^4

which is

(x^4)/5 * (pochhammer(-7,4)*C(1) + pochhammer(-6,3)*C(2) + pochhammer(-5,2)*C(3) + pochhammer(-4,1)*C(4) + pochhammer(-3,0)*C(5))

Each of the terms for x^(N-K) with order N starts with pochhammer(1-N,N-K)

In Maple notation, the expression for order N would be

`+`(seq('(`+`(seq(pochhammer(-N+K, n-K)*C[K], K = 1 .. n)))*x^(N+1-n)/(N+1-n)', n = 1 .. N))

For the general differential solution with no boundary conditions, one would add to this -c1*exp(-x) + c2 where c1 and c2 are the constants of integration.

The pochhammer(-N+K, n-K) portion can be rewritten as

GAMMA(n-N)/GAMMA(-N+K) or as (n-(N+1))!/(K-(N+1))!

Answer 2

[Teja Muppirala]

Could you define a function that is the interpolation of your discrete data points F(t)?

Fi = @(ti) interp1(t,F,ti)

This is linear interpolation, but you can choose whatever kind you want.

Now Fi is a continuous function which you can use with the ODE solver.