# Problem

How to find the initial conditions x[0],y[0] and the optimal value of parameters a,b,c,d to fit the following system of ODEs?

 x'[t] ==  a*x[t] - b*x[t]*y[t]
y'[t] == -c*y[t] + d*x[t] y[t]


Solution data is:

    (* The first column is time 't', the second column is coordinate 'x', and the last
column is coordinate 'y'. *)
data = {{11, 45.79, 41.4},    {12, 53.03, 38.9},    {13, 64.05, 36.78},
{14, 75.4, 36.04},    {15, 90.36, 33.78},   {16, 107.14, 35.4},
{17, 127.79, 34.68},  {18, 150.77, 36.61},  {19, 179.65, 37.71},
{20, 211.82, 41.98},  {21, 249.91, 45.72},  {22, 291.31, 53.1},
{23, 334.95, 65.44},  {24, 380.67, 83.},    {25, 420.28, 108.74},
{26, 445.56, 150.01}, {27, 447.63, 205.61}, {28, 414.04, 281.6},
{29, 347.04, 364.56}, {30, 265.33, 440.3},  {31, 187.57, 489.68},
{32, 128., 512.95},   {33, 85.25, 510.01},  {34, 57.17, 491.06},
{35, 39.96, 462.22},  {36, 29.22, 430.15},  {37, 22.3, 396.95},
{38, 16.52, 364.87},  {39, 14.41, 333.16},  {40, 11.58, 304.97},
{41, 10.41, 277.73},  {42, 10.17, 253.16},  {43, 7.86, 229.66},
{44, 9.23, 209.53},   {45, 8.22, 190.07},   {46, 8.76, 173.58},
{47, 7.9, 156.4},     {48, 8.38, 143.05},   {49, 9.53, 130.75},
{50, 9.33, 117.49},   {51, 9.72, 108.16},   {52, 10.55, 98.08},
{53, 13.05, 88.91},   {54, 13.58, 82.28},   {55, 16.31, 75.42},
{56, 17.75, 69.58},   {57, 20.11, 62.58},   {58, 23.98, 59.22},
{59, 28.51, 54.91},   {60, 31.61, 49.79},   {61, 37.13, 45.94},
{62, 45.06, 43.41},   {63, 53.4, 41.3},     {64, 62.39, 40.28},
{65, 72.89, 37.71},   {66, 86.92, 36.58},   {67, 103.32, 36.98},
{68, 121.7, 36.65},   {69, 144.86, 37.87},  {70, 171.92, 39.63},
{71, 202.51, 42.97},  {72, 237.69, 46.95},  {73, 276.77, 54.93},
{74, 319.76, 64.61},  {75, 362.05, 81.28},  {76, 400.11, 105.5},
{77, 427.79, 143.03}, {78, 434.56, 192.45}, {79, 410.31, 260.84},
{80, 354.18, 339.39}, {81, 278.49, 413.79}, {82, 203.72, 466.94},
{83, 141.06, 494.72}, {84, 95.08, 499.37},  {85, 66.76, 484.58},
{86, 45.41, 460.63},  {87, 33.13, 429.79},  {88, 25.89, 398.77},
{89, 20.51, 366.49},  {90, 17.11, 336.56},  {91, 12.69, 306.39},
{92, 11.76, 279.53},  {93, 11.22, 254.95},  {94, 10.29, 233.5},
{95, 8.82, 212.74},   {96, 9.51, 193.61},   {97, 8.69, 175.01},
{98, 9.53, 160.59},   {99, 8.68, 146.12},   {100, 10.82, 131.85}};


I have tried ParametricNDSolveValue but get a lot error message. Can you help me?

## Solution(up to now)

We now have a direct method to solve this problem. The main idea is interpolating many points from t=11 to t=100, then get the value of dx/dt and dy/dt by Finite Difference method(with minor time step, this can decrease the loss of accuracy). Finally, just use FindFit twice to get parameters a, b and c, d. So that system of ODEs can be determined, but how can we find x[0] and y[0]? Just solve the system backward using NDSolve or NDSolveValue.

{t, x, y} = Transpose[data]; (* extract data *)
{xdata, ydata} = (Transpose[{t, #}]) & /@ {x, y}; (* get {t,x} and {t,y} pairs *)
fx = Interpolation[xdata, Method -> "Spline"];  (* spline interpolation *)
fy = Interpolation[ydata, Method -> "Spline"];
dx=fx'[t];
dy=fy'[t];
ab=FindFit[Transpose[{fx[t], fy[t], dx}], a*X - b*X Y, {a, b}, {X, Y}]
cd=FindFit[Transpose[{fx[t], fy[t], dy}], -c*Y + d*X Y, {c, d}, {X, Y}]


Here we get {a -> 0.213493, b -> 0.00119763} and {c -> 0.104194, d -> 0.000950553}. As I put above, use NDSolve to find x[0] and y[0].

(* search backward *)
NDSolve[{X'[u] == a*X[u] - b*X [u] Y[u] /. ab,  Y'[u] == -c*Y[u] + d*X[u] Y[u] /. cd,
X[11] == 45.79, Y[11] == 41.4}, {X[0], Y[0]}, {u, 0, 12}]


It returns {{X[0] -> 10.415, Y[0] -> 102.984}}.

Let's see the plot:

(* find solution of x[t], y[t] with NDSolveValue *)
sol = NDSolveValue[{X'[u] == a*X[u] - b*X [u]Y[u] /. ab,Y'[u] == -c*Y[u] + d*X[u]Y[u] /. cd,
X[11] == 45.79, Y[11] == 41.4}, {{u, X[u]}, {u, Y[u]}}, {u, 11, 100}]
(* plot the results *)
Show[{ListPlot[{xdata, ydata}, AxesLabel -> {"time", "value"},
PlotLegends -> {"xdata", "ydata"}, PlotStyle -> {Red, Green}],
ListLinePlot[{Table[sol[[1]], {u, 11, 100, 0.01}], Table[sol[[2]], {u, 11, 100, 0.01}]},
PlotStyle -> {Orange, Blue}, PlotLegends -> {"solution.x", "solution.y"}]}]


It is obvious that scaters at the peaks and valleys are off the curves a little. This implies us that the result is not good enough and there exists improved space of the direct method. Can you image how difficult it would be using FDM when the system of ODEs is of much higher order? We must find out an universal way to solve this kind of problem.

Here we show another way:

{tval, xval, yval} = Transpose[data]; (*extract data*)
{xdata, ydata} = (Transpose[{t, #}]) & /@ {xval, yval }; (*get {t,x} and {t,y} pairs*)
(* define the embeded function with arguments a,b,c,d *)
px[a_?NumberQ, b_?NumberQ, c_?NumberQ, d_?NumberQ] := (modle[a, b, c, d] =
First[x /. NDSolve[{x'[t] == a x[t] - b x[t] y[t], y'[t] == -c y[t] + d x [t] y[t],
x[11] == 45.79, y[11] == 41.70}, x, {t, 10, 100}]]);
FindFit[xdata, px[a, b, c, d][t], {{a, .1}, {b, .1}, {c, .1}, {d, .1}}, t]
py[a_?NumberQ, b_?NumberQ, c_?NumberQ, d_?NumberQ] := (model[a, b, c, d] =
First[x /. NDSolve[{x'[t] == a x[t] - b x[t] y[t], y'[t] == -c y[t] + d x [t] y[t],
x[11] == 45.79, y[11] == 41.70}, y, {t, 10, 100}]]);
FindFit[ydata, py[a, b, c, d][t], {{a, .1}, {b, .1}, {c, .1}, {d, .1}}, t]


Unfortulately, it doesn't work. We need a good initial value for a, b, c, d, or will fail to do the fitting work. Maybe we could take a chance to get good value with the help of NSolve or FindInstance.

-
This might be stupid, but - off the top of my head... Instead of solving the ODEs system, since we have approximated values of x[t] and y[t], from which approximated values of the derivatives can be computed, why not substituting these values four at the time to solve for the four unknowns a,b,c and d? Will the average solutions found over the whole sweep of values be a sensible guess for a, b, c and d? –  Peltio Oct 25 '13 at 12:45
@Peltio I like that, even if the guess isn't good enough it will be a great starting value for FindRoot –  ssch Oct 25 '13 at 13:08
@Peltio This problem is just a simple demo for solving parameters fitting of ODEs with ParametricNDSolveValue. With data set above, we can get parameters a,b,c,d by finding a analytic solution(if exists) or simply compute derivatives by difference for substitution into those equations above. –  Kimist Oct 25 '13 at 13:27
I am not familiar with ParametricNDSolveValue, but this kind of problems can suffer from lack of precision in the data. After all predator-prey systems are used to exemplify the sensitivity of solutions from initial data, right? With my silly method I get a wide range of values for the parameters, and my tentative solution (a=0.222, b=0.01, c=0.1, d=0.001) is probably way off. –  Peltio Oct 25 '13 at 13:56
Well, my method is not that silly after all. I refined it to get {a -> 0.214, b -> 0.00122, c -> 0.105, d -> 0.00095}. If you solve the system with this parameters and plot it along with the points you have given you will find that on 0-100 it is a very good fit. Initial values are not integers though (x[0]=10.55 and y[0]=103.84), but hey, there could always be a lion with two heads and a gazelle with three legs. –  Peltio Oct 25 '13 at 18:20

This is a 'brute force' method but it seems to work. I start with extracting the data:

{tvals, xvals, yvals} = Transpose[data];
{xdata,ydata} = (Transpose[{tvals, #}])& /@ {xvals, yvals};


I then compute their interpolating functions, since I will be needing an approximate value for the derivatives:

xf = Interpolation[xdata];
yf = Interpolation[ydata];


Evaluating this will show that the derivative of the interpolations are somewhat rough

Plot[{xf'[t], yf'[t]}, {t, 11, 100}, Frame -> True, PlotRange -> All]


Since we know that this sort of problems have smooth solutions, a possible improvement could be to change xf' and yf' with a smoothed version that might be affected by less error, but I won't do that now.

The brute force idea consists in generating sets of four linear equations in the coefficients a, b, c and d and then - if the system admits a unique solution - get these solutions to get an idea about their actual values. Evaluating the system of ODEs

coeffSystem[t_] := {
xf'[t] == a*xf[t] - b*xf[t]*yf[t],
yf'[t] == -c*yf[t] + d*xf[t] yf[t]
};


at two subsequent times t1 and t2 will produce the four linear algebraic equations we need.

solvableSys[t1_, t2_] := Flatten[coeffSystem /@ {t1, t2}]


I have named this procedure solvableSys because I intended to embed a check on the determinant to make sure the system would produce a unique solution but then, it's faster to solve anyway and then discard what we do not like. A possible improvement would be selecting only 'good' systems if this can be done in a way to speed thing up.

Now, the brute force I am using here is probably an overkill but it help in producing good looking graphs for the coefficient values :-) I am solving for all pairs of t1, t2. This will produce a very long list of substitution rules

paramRules = Table[
Flatten[Solve[solvableSys[t1, t2]]], {t1, 11, 99}, {t2, t1 + 1, 100}];


from which I will extract the values for a, b, c, and d (the last substitution is to take care of the unsolvable or undetermined systems - the rude method I use to infer the solution is not affected by a bunch of added zeros).

{rawAdata, rawBdata, rawCdata, rawDdata} =
(Flatten[# /. paramRules] /. # -> 0) & /@ {a, b, c, d};


If you plot these lists of values you will see that there is a somewhat large dispersion but the values are insisting on a given line. That's what I elected to be the actual value of the parameter. Plotting the sorted data can give a clearer view of what is happening

{ {ListPlot[ Sort[rawAdata] ], ListPlot[ Sort[rawBdata] ]},
{ListPlot[ Sort[rawCdata] ], ListPlot[ Sort[rawDdata], ]} }


At this point we can extract the estimate of our parameters by inspection just by locating the inflection points and reading the parameter value on the vertical axis (the horizontal axis has no particular meaning). Doing this will get you

paramEstimate = {a -> 0.214, b -> 0.00122, c -> 0.105, d -> 0.00095}


(A better method would be computing the interpolating functions of the sorted values, and then solving for the value where their derivative is zero. The corresponding ordinate is the solution). The solution of the system on 0-100, computed with conditions at t=11, the first available point in the data, will give:

Clear[x, y];
{x[t_], y[t_]} = {x[t], y[t]} /. NDSolve[{
x'[t] == a*x[t] - b*x[t]*y[t],
y'[t] == -c*y[t] + d*x[t]*y[t],
x[11] == 45.79, y[11] == 41.4 } /. paramEstimate,
{x[t], y[t]}, {t, 0, 100}][[1]];


This solution appears to be a good fit of the supplied data over the interval 0-100.

pdata = ListPlot[{xdata, ydata}];
solved = Plot[{x[t], y[t]}, {t, 0, 100}];
Show[pdata,solved]


The initial values at t=0 are

{x[0], y[0]}
(* 10.555, 103.839 *)

-
tt = data[[All, 1]];
x = data[[All, 2]];
y = data[[All, 3]];
(* finite difference approximation of derivative *)
dx = (RotateLeft[#1] - #1 & )[x];
dy = (RotateLeft[#1] - #1 & )[y];
(* Use the least-squares optimization *)
err[a_, b_, c_, d_] = Sum[
(dx[[t]] - a*x[[t]] + b*x[[t]]*y[[t]])^2 +
(dy[[t]] + c*x[[t]] - d*x[[t]]*y[[t]])^2,
{t, -1 + Length[tt]}];
(* Set the gradient (with respect to parameters) to zero and solve *)
NSolve[{D[err[a, b, c, d], a] == 0, D[err[a, b, c, d], b] == 0,
D[err[a, b, c, d], c] == 0, D[err[a, b, c, d], d] == 0}, {a, b, c, d}]
(* Output *)
{{a -> 0.197193, b -> 0.00117457, c -> -0.0562289, d -> 0.000277672}}


I just saw that you also wanted x0 and y0.

{X[t_], Y[t_]} = {X[t], Y[t]} /.
NDSolve[{X'[t] == a*X[t] - b*X[t]*Y[t],
Y'[t] == -c*Y[t] + d*X[t]*Y[t],
X[11] == 45.79,
Y[11] == 41.4} /. {{a -> 0.197193, b -> 0.00117457,c -> -0.0562289, d -> 0.000277672}}
,{X[t], Y[t]}, {t, 0, 11}][[1]];

{X[0],Y[0]}
(* {7.67928, 20.8671} *)

-
The sum doesn't include the last element in the list because the forward finite difference is not defined there. Last@dx==Last@x-First@x is meaningless –  Timothy Wofford Oct 25 '13 at 13:56
This is just an aside that has nothing to do with the problem at hand but... your first three lines of code can be compacted with {tt,x,y}=Transpose[data]; –  Peltio Oct 25 '13 at 13:59
@Peltio, good point. I usually only think of putting one symbol on the lhs... –  Timothy Wofford Oct 25 '13 at 14:09
Correction. Last@dx==First@x-Last@x is meaningless. –  Timothy Wofford Oct 25 '13 at 14:13
@TimothyWofford Good job. But I think it is somewhat rough to use FDM(finite difference method) because of big time step. You method(least-squares optimization) can be improved by findfit. –  Kimist Oct 26 '13 at 5:50