Find parameters of ODEs to fit solution data

Question

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.

Answer 1

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 *)

Answer 2

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} *)