If you try to evaluate your model for numerical parameters, you'll notice it doesn't evaluate at all. This is because you have a typo in the definition of your model: there is one comma too many in between two of the k3r
and k19
parameters. Once you fix that, you'll also notice that the model returns a list instead of a number, which is another problem you should fix (probably with First @ NDSolve[...]
).
Finally, the way you call Model[187/100000, Amp, k1, k1r, k3, k3r, k19, k20][t]
in FindFit
will not work, since Model[...]
will evaluate to an expression involving t
and you then glue another t
at the end. Instead, Model[...]
should be returning a Function
that accepts t
as an argument. You can do that by replacing:
(Amp (InitialConditions - 1/2 (\[Beta][t] + \[Gamma][t]) + 3/2 \[Delta][t])
with:
Function[t, (Amp (InitialConditions - 1/2 (\[Beta][t] + \[Gamma][t]) + 3/2 \[Delta][t])]
and changing the solution variables in NDSolve
from {\[Beta][t], \[Gamma][t], \[Delta][t], \[Zeta][t], \[Eta]}
to {\[Beta], \[Gamma], \[Delta], \[Zeta], \[Eta]}
. Doing so will make NDSolve
return functions instead of expressions in t
.
A few other remarks: it's usually better to use NumericQ
for these kinds of matching patterns instead of NumberQ
, since NumberQ
has a few quirks that are often undesirable. For example, NumberQ[Pi]
is False
. Also, for fitting a function, you may want ask yourself if your model function should be using memoization, since FindFit
shouldn't be asking for the same value twice anyway.
The whole model will now look like:
ClearAll[Model]
Model[InitialConditions_?NumericQ, Amp_?NumericQ, k1_?NumericQ,
k1r_?NumericQ, k3_?NumericQ, k3r_?NumericQ, k19_?NumericQ,
k20_?NumericQ] :=
Function[t, Amp (InitialConditions - 1/2 (\[Beta][t] + \[Gamma][t]) + 3/2 \[Delta][t])] /. First@NDSolve[{
\[Beta]'[t] == -(k1r \[Gamma][t] + k3 (InitialConditions - 1/2 (\[Beta][t] + \[Gamma][t]) +
3/2 \[Delta][t])) \[Beta][t] + k1 (InitialConditions - 1/2 (\[Beta][t] + \[Gamma][t]) +
3/2 \[Delta][t]) + k3r \[Delta][t],
\[Gamma]'[t] == k1 (InitialConditions - 1/2 (\[Beta][t] + \[Gamma][t]) +
3/2 \[Delta][t]) - (k1r \[Beta][t] + k19 ( M - \[Zeta][t] - \[Eta][t]) +
k20 \[Zeta][t]) \[Gamma][t], \[Delta]'[t] == -k3r \[Delta][t] +
k3 \[Beta][ t] (InitialConditions - 1/2 (\[Beta][t] + \[Gamma][t]) + 3/2 \[Delta][t]),
\[Zeta]'[t] == k19 ( M - \[Zeta][t] - \[Eta][t]) \[Gamma][t] - k20 \[Zeta][t] \[Gamma][t],
\[Eta]'[t] == k20 \[Zeta][t] \[Gamma][t],
\[Beta][0] == Sqrt[(k1 InitialConditions)/k1r], \[Gamma][0] == Sqrt[(k1 InitialConditions)/k1r], \[Delta][0] == 0, \[Epsilon][0] == M, \[Zeta][0] == 0,
\[Eta][0] == 0},
{\[Beta], \[Gamma], \[Delta], \[Zeta], \[Eta]}, {t, 0, runtime}];
Test that it works:
Model[187/100000, Amp, k1, k1r, k3, k3r, k19, k20] /. Rule @@@ {{Amp, 424.5993497}, {k1, 82}, {k1r, 3000000}, {k3, 8.8*10^7}, {k3r, 1.5*10^9}, {k19, 114000}, {k20, 250}}
%[1.]
This should return a function and a number.