Given motion in response to random forcing, can a "simple" function be used to estimate future motion. Starting with the damped harmonic oscillator, we can write
$${d^2 \over {dt^2}} h(t) + 2 \zeta \omega_n {d \over dt} h(t) + \omega_n^2 h(t) = { F(t) \over m} $$
This can be rewritten in state transition form as $$ \begin{bmatrix} { d \over dt} h(t) \\ h(t) \end{bmatrix} = {e^{-\zeta \omega_n t} \over \omega_d} \begin{bmatrix} \omega_d cos(\omega_d t) + \zeta \omega_n sin(\omega_d t) & sin(\omega_d t) \\ -\omega_n^2 sin(\omega_d \tau) & \omega_d cos(\omega_d t) - \zeta \omega_n sin(\omega_d t) \end{bmatrix} \begin{bmatrix} { d \over dt} h(0) \\ h(0) \end{bmatrix} \\ + {1 \over {m \omega_d}} \int_0^\tau e^{-\zeta \omega_n t} \begin{bmatrix} -sin(\omega_d \tau) \\ - \omega_d cos(\omega_d \tau) + \zeta \omega_n sin(\omega_d \tau) \end{bmatrix} F(t-\tau) d \tau $$
How to approximate the integral term. The goal is to use previous motion data to estimate the parameters $\omega$ and $\zeta$ and generate an approximation for F and then forecast future data. It may be simpler to assume that F was harmonic. The convolution integral can probably be used to split the integral into two manageable parts, which is what I was trying to do, before posting this question. An alternate approach would be to assume F is a sum of a small number of harmonics, and fit F to the data.
If we assume the forcing is of the form $$C_i sin(Omega_i t + phi_i)$$ the a closed form solution can be found.
mupad commands
assume(z>0.001 and z<.3), assume(omega>0 and omega<.25)
IVP := ode({h''(t) +2*z*omega*h'(t)+ omega*omega*h(t)=1/m*C_i*sin(Omega_i*t+phi_i), h(0) = h0, h'(0) = dh0}, h(t))
sol := solve(IVP)
I'm not sure how to simplify the solution further, but this solution could be used to estimate the future values of h from initial values of h and dh.