How to fit and ODE to data?

Question

Consider the following ODE $$ y'(t)=\alpha x(t)-\beta y(t) $$ and the following datasets $$ X=\{(t_0,x_0),...,(t_n,x_n)\}\\ Y=\{(t_0,y_0),...,(t_n,y_n)\} $$ How can I find $\alpha$ and $\beta$ that best fits this data? In particular, is it possible to use the following solution $$ y(t) = e^{-\beta t} \left( \alpha \int_0^t e^{\beta \tau} x(\tau) \, d\tau + C \right), $$ directly, where $C$ is a constant?

Thoughts: This a fairly new field for me and I have briefly heard about inverse problems in ODEs, but rather than setting up a global minimizer that numerically finds it, I am wondering whether there is an analytical way of solving this problem, to the best possible fit in some sense.

Example: As an example, one can consider the following data $$ \begin{align} X&=\{91, 110, 125, 105, 88, 84\}\\ Y&=\{1.0, 0.97, 1.0, 0.95, 0.92, 0.89\} \end{align} $$ at time points $$ T=\{0, 5, 9, 18, 28, 38\} $$ Using a numerical approach (see Python script here), I get the following fit

with $\alpha\simeq 0.000486$, $\beta\simeq 0.057717$, and a chi-square error of $0.0017296$. Is it possible to improve this? I have tried global optimizers, but to no avail, hence thinking of a potential analytical approach.

Answer 1

There is too little data to be able to come up with a decent fit. With the little data available, a fit must be made with few parameters to determine. In this sense, we first interpolate the data to $x(t)$ in spline form and continue with obtaining $\alpha$, $\beta$ and $y_0$ with a nonlinear least squares fit, on the solution of the integral that in this case can be obtained explicitly. The ideal would be to carry out the adjustment of $x(t)$ and $y(t)$ simultaneously, but for this, a more significant data set would be necessary. Follows a MATHEMATICA script to accomplish that.

tk = {0, 5, 9, 18, 28, 38};
xk = {91, 110, 125, 105, 88, 84};
yk = {1, 0.97, 1, 0.95, 0.92, 0.89};
base[t_] := {1, t, t^2}
n = Length[tk];
X[t_] := Sum[base[t - tk[[k]]] . {a0[k], a1[k], a2[k]} (UnitStep[t - tk[[k]]] - UnitStep[t - tk[[k + 1]]]), {k, 1, n - 1}]
equsa1 = Table[base[tk[[k + 1]] - tk[[k]]] . {a0[k], a1[k], a2[k]} - a0[k + 1], {k, 1, n - 1}];
subsa1 = Solve[equsa1 == 0, Table[a2[k], {k, 1, n - 1}]][[1]];
equsa2 = Table[a1[k] + 2 a2[k] (tk[[k + 1]] - tk[[k]]) - a1[k + 1], {k, 1, n - 1}];
subsa2 = Solve[equsa2 == 0, Table[a2[k], {k, 1, n - 1}]][[1]];
equstot = Join[Join[equsa1, equsa2]];
varstot = Join[Join[Table[a1[k], {k, 1, n - 1}], Table[a2[k], {k, 1, n - 1}]]];
sols = Solve[equstot == 0, varstot][[1]];
Xt = X[t] /. sols;
subsvars = Table[a0[k] -> xk[[k]], {k, 1, n}];
Xt0 = Xt /. subsvars /. {a1[n] -> 0};
gr1 = ListPlot[Transpose[{tk, xk}], PlotStyle -> {Red, PointSize[0.02]}];
Show[gr1, Plot[Xt0, {t, tk[[1]], tk[[n]]}], PlotRange -> All]
gr2 = ListPlot[Transpose[{tk, yk}], PlotStyle -> {Red, PointSize[0.02]}];

tmax = tk[[n]]
Clear[x]
soly = DSolve[{y'[t] == alpha Xt0 - beta y[t], y[0] == y0}, y, {t, 0, tmax}][[1]];
yt = y[t] /. soly;
obj = Sum[((yt /. {t -> tk[[k]]}) - y2[[k]])^2, {k, 1, n}];
sol = Minimize[obj, {alpha, beta, y0}]

(** {0.000827948, {alpha -> 0.0000672746, beta -> 0.0102676, y0 -> 1.00217}} **)

yt0 = yt /. sol[[2]];
gr3 = Plot[yt0, {t, 0, tmax}];
Show[gr3, gr2]

Answer 2

They are two data sets $(x,t)$ and $(y,t)$ to be fitted to two functions $x(t)$ and $y(t)$ respectively. Thus we need two equations in which some parameters have to be optimised.

In the wording of the problem only one equation is specified : $$y'(t)=\alpha x(t)-\beta y(t)\tag1$$ in which two parameters $\alpha$ and $\beta$ have to be optimised.

So one equation is missing. That is why (for example) a polynomial equation might be proposed to approximate $y(t)$. With this additional equation the problem could be solved thanks to numerical calculus. But choosing a polynomial is rather arbitrary and the result will depend on the degree.

A different approach is proposed. This doesn't entirely solve the problem but leads directly to rough approximates of $\alpha$ and $\beta$.

Integrating the given equation $(1)$ leads to :

$$y(t)=y_0+\alpha \int_{t_0}^t x(\tau)d\tau-\beta \int_{t_0}^t y(\tau)d\tau\tag2$$ Approximates of the integrals can be computed thanks to numerical integration :

$$S_0=0\quad;\quad S_i=S_{i-1}+\frac12(x_i+x_{i-1})(t_i-t_{i-1})\quad \text{ from } i=1 \text{ to } i=n.$$ $$S_i\simeq\int_{t_0}^{t_i} x(\tau)d\tau$$

$$R_0=0\quad;\quad R_i=R_{i-1}+\frac12(y_i+y_{i-1})(t_i-t_{i-1})\quad \text{ from } i=1 \text{ to } i=n.$$ $$R_i\simeq\int_{t_0}^{t_i} y(\tau)d\tau$$ Approximate of Eq.$(2)$ gives the linear system : $$y_i-y_0\simeq \alpha S_i-\beta R_i \quad \text{ from } i=0 \text{ to } i=n.\tag3$$ Linear regression gives approximates of $\alpha$ and $\beta$.

Of course this doesn't entirely solve the problem as pointed out at the begining of the answer. The functions $x(y)$ and $y(x)$ are not explicitely found. The modeling of the experiment should be completed. If not, one have to guess an additional equation possibly a convenient kind of function for $x(t)$ including parameters to optimise. This would be the most difficult part of the task even if approximates of $\alpha$ and $\beta$ are already known.

IN ADDITION : Example of fitting with a polynomial.

The choice of polynomials to approximate $x(t)$ and $y(t)$ might be not the best but the calculus is much simpler than with a no-elementary function.

Starting with the approximate $\alpha$ and $\beta$ computed above a linear regression leads to the coefficients of the polynomial.

Note : The magnitude of $x$ and $y$ are very different. It is simpler to choose LMSRE instead of LMSE as criteria of global fitting. This avoid rescaling.

Example :

$p =$ degree of the polynomial. (For example $p=5$).

$n =$ Number of points. (For example $n=6$).