Parameter identification: fitting data to a model that is defined using NDSolve

You can, in principle, fit any data to any model by constructing an expression for the difference between the data and the model, and using  to minimize that expression with respect to the parameters of the model. The example below shows one way to do this when the model is computed by , and the difference to be minimized is the sum of the squared differences between the data and the model at each data point (least-squares fitting).

If the model is given by a function f, the sum of the squared differences between the data and the model can be constructed using

\(data\  = 
        \ {{0.5, \ 6.2}, \ {1.1, \ 8.1}, \ {1.5, \ 8.8}, 

          \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {2.1, \ 8.3}, 
          \ {2.3, \ 6.8}, \ {3.1, \ 2.9}}\ ; \)

ss=Plus @@ Apply[(f[#1]  -  #2)^2 &,  data,  {1}]

In this example the function f is the solution of a differential equation

f' (x)={({{10 - f(x),  x < p}, {-f(x),  x ≥ p}})
f(0)=q

ie.

eqn = {f'[x]  ==  If[x  <   p,  10  -  f[x],  -f[x]], 
         f[0]  ==  q}

where p and q are parameters of the solution.

A function sse[p, q] to compute the sum of the squared differences between the data and the model in this example can be written as

\(sse[p_?NumberQ, \ q_?NumberQ]\  := \ 
        Block[{sol, \ f}, 
		
          sol\  = \ \(NDSolve[eqn, \ f, \ {x, \ 0, \ 5}]\)[\([1]\)]; 
          ss\  /. \ sol]\)

This expression can be minimized using  to obtain the best fit values for p and q.

FindMinimum[sse[p,  q],  {p, 3,  5},  {q, 23, 20}, 
      MaxIterations -> 100]

These calculations can be computationally and mathematically quite demanding. In addition to all of the usual difficulties with non-linear fitting, solving a differential equation is a non-trivial task, and when the model is defined by a differential equation, that differential equation must be solved (with new parameters) at each step in the minimization process.

Another example:

Quit[]

In this example the function f is the solution of a differential equation

f' (t) = a f(t) (b - f(t))f(0) = c

ie.

eqn = {f'[t]  ==  f[t]*(b  -  a*f[t]) ,  f[0]  ==  c}

where a, b, and c are parameters of the solution.

data1 = {{1.30, 106}, {1.40, 123}, {1.50, 132}, {1.60, 151}, {1.70, 
          179}, {1.80, 205}, {1.90, 226}};

ss=Plus @@ Apply[(f[#1]  -  #2)^2 &,  data1,  {1}]

\(sse[a_?NumberQ, \ b_?NumberQ, \ c_?NumberQ]\  := \ 
        Block[{sol, \ f}, 
		
          sol\  = \ 
            \(NDSolve[{\(f'\)[t]\  == \ f[t]*\((b\  - a*f[t])\)\ , \ 
                  f[0]\  == \ c}, \ f, \ {t, 0, \ 1.90}, MaxSteps -> 2000, 
                AccuracyGoal -> 10]\)[\([1]\)];  (*\ 
            Print["a=", a, ", b=", b, ", c=", c]; *) \ 
          ss\  /. \ sol]\)

a0=0.001;b0=1.23;c0=100;
Do[
min=FindMinimum[sse[a,b,c], {a,a0,2*a0}, {b,b0,2*b0},{c,c0,2*c0},
    MaxIterations->1000,AccuracyGoal->10,PrecisionGoal->10,
    WorkingPrecision->30];
a0=a/.min[[2]];b0=b/.min[[2]];c0=c/.min[[2]]; {2}]; 
min

dplot = ListPlot[data1]

\(\(	sol\  = \ 
      \(NDSolve[{\(f'\)[t]\  == \ a*f[t]*\((b\  - \ f[t])\)\ , \ 
              f[0]\  == \ c} /. min[\([2]\)], \ f, \ {t, 0, \ 1.90}, 
          MaxSteps -> 2000, AccuracyGoal -> 10]\)[\([1]\)]\)\)

splot = Plot[f[s] /. First[sol], {s, 0, 1.9}]

Show[dplot, splot]

plt = 
      Plot3D[Evaluate[Log[sse[a, b, 19.5]]], {a, 0.00001, 0.0009}, {b, 1, 2}, 
        PlotRange -> All, PlotPoints -> 40]

Show[plt, ViewPoint -> {0,  -3.383,  0}]