MATLAB Answers

Solving Second Order Differential Equations with Discrete Time Data

33 views (last 30 days)
Kevin on 13 May 2011
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.

Answers (2)

Walter Roberson
Walter Roberson on 13 May 2011
Is that a system of equations,
g''(t) + g'(t) = F(t)
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
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
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))!

Sign in to comment.

Teja Muppirala
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.

Community Treasure Hunt

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

Start Hunting!