# Solving ODE System and Estimating Parameters with Experimental Data

I am new to Mathematica. I am trying to find out the best-fit values of the parameters appearing in an ODE system of three equations. The three ODE equations are:

{a'[t] == -k[1]*a[t], b'[t] == -k[2]*b[t] + k[1]*a[t], c'[t] == k[2]*b[t]}


The experimental data set for the three equations with respect to time (in minutes) is:

adata = {{0, 1}, {2, 0.88}, {6, .69}, {10, .53}, {20, .28}, {30, .15}, {50, .043}, {70, .012}, {90, 0}, {120, 0}, {150, 0}, {200, 0}}
bdata = {{0, 0}, {2, 0.12}, {6, .29}, {10, .42}, {20, .56}, {30, .57}, {50, .46}, {70, .33}, {90, 0.22}, {120, 0.12}, {150, .06}, {200, 0.02}}
cdata = {{0, 0}, {2, 0.003}, {6, .030}, {10, .050}, {20, .16}, {30, .28}, {50, .50}, {70, .66}, {90, 0.78}, {120, 0.88}, {150, .94}, {200, 0.98}}


(First entry in each pair is time)
The goal here is to solve ODE, and then optimize the parameters (k1, k2) such that the equations can best describe the experimental data. I would be very thankful to you, if you can spare some of your valuable time and help me with this query. the objective fuction to be minimzed is as follows

But it is giving me the error like 'not a real number' etc.

• Welcome to Mathematica.SE! I have edited your post the way we like to see things. Note that we have entered all equations and data in properly formatted code blocks with proper Mathematica syntax. This way, we can copy and paste your code into our own copies of Mathematica without having to write it all out ourselves. Feb 19, 2020 at 16:47
• @march thank you for editing my post. Feb 20, 2020 at 9:16
• @morbo, or any body else could you please see why the Functionn 'SS' is giving error when evaluated through FindMinimum or NMinimize? Feb 20, 2020 at 14:55

ClearAll["Global*"]

adata = {{0, 1}, {2, 0.88}, {6, .69}, {10, .53}, {20, .28}, {30, .15}, {50, .043}, {70, .012}, {90, 0}, {120, 0}, {150, 0}, {200, 0}};
bdata = {{0, 0}, {2, 0.12}, {6, .29}, {10, .42}, {20, .56}, {30, .57}, {50, .46}, {70, .33}, {90, 0.22}, {120, 0.12}, {150, .06}, {200, 0.02}};
cdata = {{0, 0}, {2, 0.003}, {6, .030}, {10, .050}, {20, .16}, {30, .28}, {50, .50}, {70, .66}, {90, 0.78}, {120, 0.88}, {150, .94}, {200, 0.98}};

eqns = {a'[t] == -k1 a[t], b'[t] == -k2 b[t] + k1 a[t],  c'[t] == k2 b[t]};
Thread[{aa[t_], bb[t_], cc[t_]} = DSolveValue[{eqns, a[0] == 1, b[0] == 0, c[0] == 0}, {a[t], b[t],  c[t]}, t]];

model[k1_, k2_] := Sum[(aa[adata[[i, 1]]] - adata[[i, 2]])^2 + (bb[bdata[[i, 1]]] -
bdata[[i, 2]])^2 + (cc[cdata[[i, 1]]] - cdata[[i, 2]])^2, {i, Length@adata}]

fit = Last@ NMinimize[model[k1, k2], {k1, k2}]


{k1 -> 0.0632016, k2 -> 0.0211064}

Thread[{k1, k2} = Values@fit];

Show[Plot[{aa[t], bb[t], cc[t]}, {t, 0, 200}, Frame -> True], ListPlot[{adata, bdata, cdata}]]


• Thank you @OkkesDulgerci, That was very helpful. Feb 21, 2020 at 9:35

First thing first, when you have data points, always plot and have a look at it.

I replaced your k[1] and k[2] with $$\delta$$ and $$\alpha$$.

eqns = {a'[t] == -k[1]*a[t], b'[t] == -k[2]*b[t] + k[1]*a[t], c'[t] == k[2]*b[t]} /. {k[1] -> \[Delta], k[2] -> \[Alpha]}

ListLinePlot[{adata, bdata, cdata}, ImageSize -> Large]


Looks like a bunch of $$a e^{x}$$ functions to me...And your system of equations is small, they may have an analytical solution.

DSolve[{eqns, a[0] == 1, b[0] == 0, c[0] == 0}, {a[t], c[t], b[t]}, t]


$$\left\{\left\{a(t)\to e^{\delta (-t)},b(t)\to -\frac{\delta e^{\alpha (-t)-\delta t} \left(e^{\delta t}-e^{\alpha t}\right)}{\alpha -\delta },c(t)\to \frac{e^{\alpha (-t)-\delta t} \left(\alpha e^{\alpha t+\delta t}-\delta e^{\alpha t+\delta t}-\alpha e^{\alpha t}+\delta e^{\delta t}\right)}{\alpha -\delta }\right\}\right\}.$$

Oh, they do. We're in luck! We can write a simple manipulate for fitting.

Manipulate[Plot[{E^(-t \[Delta]), -((E^(-t \[Alpha] - t \[Delta]) (-E^(t \[Alpha]) + E^(t \[Delta])) \[Delta])/(\[Alpha] - \[Delta])), (E^(-t \[Alpha] - t \[Delta]) (-E^(t \[Alpha]) \[Alpha] + E^(t \[Alpha] + t \[Delta]) \[Alpha] + E^(t \[Delta]) \[Delta] - E^(t \[Alpha] + t \[Delta]) \[Delta]))/(\[Alpha] - \[Delta])}, {t, 0.1, 199}, Epilog -> {{Red, Point[adata]}, {Blue, Point[bdata]}, {Green, Point[cdata]}}, ImageSize -> Large], {{\[Delta], 0.0613}, 0.0001, 0.08}, {{\[Alpha], 0.02085}, 0.0001, 0.03} ]


We've found by hand some good initial values for $$\delta$$ and $$\alpha$$

We can use this as initial values for a nonlinear model fit.

m1 = NonlinearModelFit[
adata, {a[t]} /. sol, {{\[Delta], 0.0613}, {\[Alpha], 0.02085}}, t]
m2 = NonlinearModelFit[
bdata, {b[t]} /. sol, {{\[Delta], 0.0613}, {\[Alpha], 0.02085}}, t]
m3 = NonlinearModelFit[
cdata, {c[t]} /. sol, {{\[Delta], 0.0613}, {\[Alpha], 0.02085}}, t]


Plot the result

Plot[{m1[t], m2[t], m3[t]}, {t, 0, 200}, PlotRange -> All,
Epilog -> {{Red, Point[adata]}, {Blue, Point[bdata]}, {Green,
Point[cdata]}}]


and then get some infos about our parameters that NonlinearModelfit found.

m1["ParameterTable"]
m2["ParameterTable"]
m3["ParameterTable"]


$$\begin{array}{l|llll} \text{} & \text{Estimate} & \text{Standard Error} & \text{t-Statistic} & \text{P-Value} \\ \hline \delta & 0.063221 & 0.000225209 & 280.722 & \text{8.092528367214482\grave{ }*{}^{\wedge}-21} \\ \alpha & 0.02085 & 0. & \infty & \text{0\grave{ }\grave{ }323.6072453387798} \\ \end{array}$$

$$\begin{array}{l|llll} \text{} & \text{Estimate} & \text{Standard Error} & \text{t-Statistic} & \text{P-Value} \\ \hline \delta & 0.0629746 & 0.000365423 & 172.333 & \text{1.0634929126057997\grave{ }*{}^{\wedge}-18} \\ \alpha & 0.0210366 & 0.0000895516 & 234.91 & \text{4.805391875551799\grave{ }*{}^{\wedge}-20} \\ \end{array}$$

$$\begin{array}{l|llll} \text{} & \text{Estimate} & \text{Standard Error} & \text{t-Statistic} & \text{P-Value} \\ \hline \delta & 0.0640477 & 0.00221151 & 28.961 & \text{5.61521578981255\grave{ }*{}^{\wedge}-11} \\ \alpha & 0.0210622 & 0.000314458 & 66.9795 & \text{1.3404461637525215\grave{ }*{}^{\wedge}-14} \\ \end{array}$$

Histogram[m1["FitResiduals"]]
Histogram[m2["FitResiduals"]]
Histogram[m3["FitResiduals"]]


• Thank you @morbo. Can you please clarify to me that this 'Manipulate' does the same job as 'FindMinimum' or 'NMinimize' in terms of regressing to get best-fit parameters? Also how do I upvote, I cant see a button? Feb 20, 2020 at 9:23
• what are the values we got for parameters δ and α? I actually need these values for further input in another software. Feb 20, 2020 at 9:31
• If you click on the plus beaide the slider you‘ll get the values i found playing around. I will answer the rest later when i‘m on my computer...$\delta=0.0631$ and $\alpha=0.02085$ Feb 20, 2020 at 9:36
• @wahab the blue arrows up and down beside my answer and the grey check mark is there too Feb 20, 2020 at 9:37
• Thank you @morbo. Done. Feb 20, 2020 at 9:54

Using the language resources

adata = {{0, 1}, {2, 0.88}, {6, .69}, {10, .53}, {20, .28}, {30, .15}, {50, .043}, {70, .012}, {90, 0}, {120, 0}, {150, 0}, {200, 0}};
bdata = {{0, 0}, {2, 0.12}, {6, .29}, {10, .42}, {20, .56}, {30, .57}, {50, .46}, {70, .33}, {90, 0.22}, {120, 0.12}, {150, .06}, {200, 0.02}};
cdata = {{0, 0}, {2, 0.003}, {6, .030}, {10, .050}, {20, .16}, {30, .28}, {50, .50}, {70, .66}, {90, 0.78}, {120, 0.88}, {150, .94}, {200, 0.98}};
bdata0 = Table[Flatten[{2, bdata[[k]]}], {k, 1, Length[adata]}];
cdata0 = Table[Flatten[{3, cdata[[k]]}], {k, 1, Length[adata]}];

tmax = 200;
eqns = {a'[t] == -k1*a[t], a[0] == a0, b'[t] == -k2*b[t] + k1*a[t], b[0] == b0, c'[t] == k2*b[t], c[0] == c0};
sol = ParametricNDSolveValue[eqns, {a, b, c}, {t, 0, tmax}, {k1, k2, a0, b0, c0}];

model[k1_, k2_, a0_, b0_, c0_][i_, t_] := Through[sol[k1, k2, a0, b0, c0][t], List][[i]] /; And @@ NumericQ /@ {k1, k2, a0, b0, c0, i, t};
fit = NonlinearModelFit[transformedData, model[k1, k2, a0, b0, c0][i, t], {k1, k2, a0, b0, c0}, {i, t}] // Quiet;

fit["ParameterTable"]


obtaining the results

$$\left[ \begin{array}{ccccc} \text{} & \text{Estimate} & \text{Standard Error} & \text{t-Statistic} & \text{P-Value} \\ \text{k1} & 0.0630922 & 0.000256218 & 246.244 & \text{1.358897\grave{ }*{}^{\wedge}-52} \\ \text{k2} & 0.0210129 & 0.0000735043 & 285.874 & \text{1.3335103\grave{ }*{}^{\wedge}-54} \\ \text{a0} & 0.99959 & 0.00160667 & 622.151 & \text{4.545291\grave{ }*{}^{\wedge}-65} \\ \text{b0} & 0.000233603 & 0.00172373 & 0.135522 & 0.893075 \\ \text{c0} & 0.00247035 & 0.00120983 & 2.04189 & 0.0497497 \\ \end{array} \right]$$

• You can use a[0]==1,b[0]==0,c[0]==0` from data.. Feb 20, 2020 at 12:33
• I know that but I prefer to consider the model as stated because it is more general and can be applied to other dynamic systems. The data could contain noise as well. Feb 20, 2020 at 12:36