2
$\begingroup$

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

enter image description here But it is giving me the error like 'not a real number' etc.

|improve this question|||||
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – march Feb 19 at 16:47
  • $\begingroup$ @march thank you for editing my post. $\endgroup$ – wahab Feb 20 at 9:16
  • $\begingroup$ @morbo, or any body else could you please see why the Functionn 'SS' is giving error when evaluated through FindMinimum or NMinimize? $\endgroup$ – wahab Feb 20 at 14:55
1
$\begingroup$ 
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}]]

enter image description here

|improve this answer|||||
$\endgroup$
  • $\begingroup$ Thank you @OkkesDulgerci, That was very helpful. $\endgroup$ – wahab Feb 21 at 9:35
2
$\begingroup$

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]

plot

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} ]

manupulate

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]}}]

new plot

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"]]

histo

|improve this answer|||||
$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – wahab Feb 20 at 9:23
  • $\begingroup$ what are the values we got for parameters δ and α? I actually need these values for further input in another software. $\endgroup$ – wahab Feb 20 at 9:31
  • $\begingroup$ 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$ $\endgroup$ – morbo Feb 20 at 9:36
  • $\begingroup$ @wahab the blue arrows up and down beside my answer and the grey check mark is there too $\endgroup$ – morbo Feb 20 at 9:37
  • $\begingroup$ Thank you @morbo. Done. $\endgroup$ – wahab Feb 20 at 9:54
0
$\begingroup$

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}};
adata0 = Table[Flatten[{1, adata[[k]]}], {k, 1, Length[adata]}];
bdata0 = Table[Flatten[{2, bdata[[k]]}], {k, 1, Length[adata]}];
cdata0 = Table[Flatten[{3, cdata[[k]]}], {k, 1, Length[adata]}];
transformedData = Join[Join[adata0, bdata0], cdata0];

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] $$

|improve this answer|||||
$\endgroup$
  • $\begingroup$ You can use a[0]==1,b[0]==0,c[0]==0 from data.. $\endgroup$ – OkkesDulgerci Feb 20 at 12:33
  • $\begingroup$ 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. $\endgroup$ – Cesareo Feb 20 at 12:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.