Solving Second Order Differential Equations with Discrete Time Data
4 次查看(过去 30 天)
显示 更早的评论
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
0 个评论
回答(2 个)
Walter Roberson
2011-5-13
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 ?
3 个评论
Walter Roberson
2011-5-13
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
2011-5-14
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
2011-5-14
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.
0 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!