# Algorithm for finding initial conditions of differential equations given trajectory

Let's say I'm given a system of three first-order differential equations in three variables, where all of the equations are known, and we additionally know the trajectory of two of the variables at a set of timepoints. Is there a fast algorithm to determine the starting value of the third variable, whose trajectory is unknown.

The best I've been able to come up with so far randomly initializing the starting point of the third variable, and then using a numerical solver to determine the hypothetical trajectory that the two variables with known trajectory follow, and then optimizing on the sum of the difference of squares between the hypothetical trajectory and the actual trajectories of these two variables. However, this can take a while to converge, which is why I'm looking for something faster.

• Your second paragraph describes a standard approach for solving this sort of problem (called 4DVAR in numerical weather prediction, where finding initial conditions from observations of the state are the crucial step in getting reasonably accurate forecasts). A very popular alternative approach is the (Ensemble) Kalman Filter. But be aware that this problem is difficult (in fact, ill-posed in general) and there is no free lunch, so don't expect to find a fast, stable, and accurate method (pick two). – Christian Clason Dec 1 '17 at 9:49
• But depending on how you optimize currently, there may be room for improvement -- for example, this problem can be solved with Newton methods. – Christian Clason Dec 1 '17 at 12:33

As Christian Clason alluded to in his comment, you have reinvented the ensemble Kalman filter! The Kalman filter is an algorithm for estimating, from imperfect measurements, the true state of a system governed by a set of first-order ODEs or difference equations. For example, you might only be able to make noisy measurements of the state of the system, or you might only be able to measure some of the state variables but not all. The measurements $y$ are assumed to be a linear function of the state variables $x$, plus some noise $\eta$:

$y = Hx + \eta$

where $H$ is the measurement matrix. In your case, the measurement matrix is not of full rank.

One of the great parts about the Kalman filter is that it gives you not just an estimate of the state of the system but a probability distribution, assuming that the errors are normally distributed, even when the measurement matrix is not of full rank. While the backward ODE may be unstable, by also tracking covariance of the Kalman filter estimates of the system state, you can also put some bound on the uncertainty of your guess for the initial condition. This doesn't remedy your lack of information; it only tells you how this lack of information, i.e. only being able to observe the first two variables, results in uncertainty in the third variable.

Kaipio and Somersalo have a chapter on Kalman filtering. I also found these notes to be really helpful. If your ODE system is nonlinear, it might be beneficial to try implementing the algorithm for linear systems of ODE first. It's especially valuable to see how the results change when the governing system is stable, unstable, or oscillatory.

Suppose your differential equations are: $$\dot{x} = f_x(x,y,z)\\ \dot{y} = f_y(x,y,z)\\ \dot{z} = f_z(x,y,z),$$ and you know $x$ and $y$ at (ordered) times $t_1, t_2, …$.

If $t_{i+1}-t_i$ is sufficiently small for some $i$ (preferably all of them), you have two ways to estimate the derivative of the first variable at : \begin{alignat*}{1} \dot{x}\left( \frac{t_i+t_{i+1}}{2}\right ) &\approx \frac{x(t_{i+1})-x(t_i)}{t_{i+1}-t_i} \\ \dot{x}\left( \frac{t_i+t_{i+1}}{2}\right ) &\approx f_x \left( \frac{x(t_i)+x(t_{i+1})}{2}, \frac{y(t_i)+y(t_{i+1})}{2}, z\left( \frac{t_i+t_{i+1}}{2}\right ) \right). \end{alignat*} If you equate the right-hand sides, the only unknown in this equation isz\left( \frac{t_{i+1}+t_i}{2}\right )$, so all you have to do is to solve it for this term. You can then apply the same approach to$y$to get a second estimate and refine your result. (While I used the middle of the interval, you can apply the same approach to the beginning or end as well. This has less accuracy but may be more comfortable and even faster, as it saves one step in the following.) Once you have this estimate for$z\left( \frac{t_{i+1}+t_i}{2}\right )$, you can integrate backwards and forwards to the known time points and refine your result. The best strategy here depends on the accuracy of your data and how your times$t_i$are distributed. If$t_1$is not your initial time, first obtain an accurate estimate of$z(t_1)$and then integrate backwards to your initial time. This may not be much different from what you have been doing, but you should have a better initial guess. Be aware of the possibility that your backwards differential equation may be unstable (e.g., for strongly dissipative systems), in which case you may want to avoid backwards integrations as far as possible in the above (or only perform them for very short time steps), in which you are essentially back at your method. Also, if$t_1\$ is much later than your initial time (with respect to instability time scales), you have the more fundamental problem that the solution to your problem is overly sensitive on measurement and numerical errors – every solution to your problem would also pose a method to backwards-integrate unstable differential equations.

• That's not a good idea in general, since the backward equation may very well be unstable... – Christian Clason Dec 1 '17 at 9:44
• the backward equation may very well be unstable – Yes, but then the asker’s problem is not really solvable anyway, irrespective of the method (also see my edit). – Wrzlprmft Dec 1 '17 at 9:53
• No, it's just more difficult (since it is ill-posed) and requires special methods. (Unless "solvable" means "solvable exactly" without requiring approximation, which is an overly strict definition for any numerical method.) But you are correct that the fundamental problem is the strong sensitivity to measurement errors. – Christian Clason Dec 1 '17 at 9:56
• @ChristianClason: Thinking about this, let me refine this (see my edit). Also, we may have run into a problem due to the several notions of stability in this context. What I have in mind right now are cases like a ball with friction being almost at rest at the bottom of a valley or a typical chaotic oscillator where you want to go backwards for a considerable number of oscillations – those are simply impossible to solve properly under experimental conditions since the possible initial states for a given final state are all over the place once you have a tiny inaccuracy. – Wrzlprmft Dec 1 '17 at 10:16
• This is exactly what I mean -- your comment has an intuitive explanation, but there is also a rigorous functional-analytic statement for this (basically, the forward operator is compact and has exponentially decaying singular values). You are also correct that making this well-posed requires imposing additional information (e.g., "if in doubt, take the lowest-energy initial state" -- again, there's a rigorous mathematical formulation of this). I'm fond of quoting Lanczos here: a lack of information cannot be remedied by any mathematical trickery. – Christian Clason Dec 1 '17 at 12:38