33 views (last 30 days)

Show older comments

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

Walter Roberson
on 13 May 2011

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 ?

Walter Roberson
on 14 May 2011

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))!

Teja Muppirala
on 14 May 2011

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.

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!