I wasn't sure how to set it up with the standard fit functions, so I just rolled my own least squares. (It probably can be done, but I had an inkling I might want to have greater control over the computation. I hope it helps.)
I start by defining m1
and m2
on your data:
Clear[m1, m2];
(m1[t_] /; t == First@# = Last@#) & /@ m1data;
(m2[t_] /; t == First@# = Last@#) & /@ m2data;
Then I defined the convolution, using the trapezoid rule and based on a consisted delta t
of 0.1
in both data sets.
ClearAll[ni];
SetAttributes[ni, Listable];
ni[0. | 0, _] := 0.;
mem : ni[t_?NumericQ, a2_?NumericQ] := mem = Quiet[
NIntegrate[m1[tau]*Exp[-a2 (t - tau)], {tau, 0, t},
Method -> {"TrapezoidalRule", "Points" -> 1 + Round[t/0.1]},
MaxRecursion -> 0],
NIntegrate::ncvb]
The objective function is the sum of squares of the residuals the solution to your ODE. This is to be minimized.
obj = m2[t] - (m1[t] + a1*ni[t, a2]) /. {m1[t] -> m1data[[All, 2]],
m2[t] -> m2data[[All, 2]], t -> m1data[[All, 1]]} // #.# &;
The value for a1
is consistent, but the model is not very sensitive to a2
. It's always a big value, which suggests only t = tau
is contributing much to the convolution.
FindMinimum[obj, {{a1, -1.}, {a2, 100}}]
(* {1.20002, {a1 -> -20.173, a2 -> 266.816}} *)
FindMinimum[obj, {{a1, -13.}, {a2, 1550}}]
(* {1.20002, {a1 -> -20.173, a2 -> 1550.}} *)
FindMinimum[obj, {{a1, -20.}, {a2, 50}}]
(* {1.20002, {a1 -> -20.173, a2 -> 2250.22}} *)
FindMinimum[obj2, {{a1, -53.}, {a2, 4500}}]
(* {1.20002, {a1 -> -20.173, a2 -> 4500.}} *)
Note: The trapezoid rule in NIntegrate
uses Romberg quadrature. If you turn it off, then a1
is consistently -13.4487
and a2
is very large just as before.
m1data
? $\endgroup$m1[tau]
orm1[t]
? If the former, do you not need to give the definition ofm1[tau]
in the solution of the ODE ? $\endgroup$