# Mathematica unable to interpret equations in DSolve

I'm trying to solve a system of seven ODEs with seven conditions.

Since the following approach using DSolve is yielding the same output as the input, Mathematica can't find a solution.

DSolve[
{EA'[t] == k1 e[t] a[t] + k22 EAB[t] - k2 b[t] EA[t] - k11 EA[t],
EB'[t] == k4 e[t] b[t] - k44  EB[t] + k33 EAB[t] - k3 a[t] eb[t],
EAB'[t] ==
k2 EA[t] b[t] + k3 EB[t] a[t] - k22 EAB[t] - k5 EAB[t] + k33 EAB[t],
p'[t] == k5 EAB[t],
e'[t] == k44 EB[t] + k11 EA[t] - k4 e[t] b[t] - k1 e[t] a[t] ,
a'[t] == k11 EA[t] - k1 e[t] a[t] + k33 EAB[t] - k3 EB[t] a[t],
b'[t] == k44 EB [t] - k4 e[t] b[t] + k22 EAB[t] - k2 EA[t] b[t] ,

EA[0] == 0,
EB[0] == 0,
EAB[0] == 0,
a[0] == 1,
e[0] == 1/10,
b[0] == 100,
p[0] == 0
}, {EAB[t], EB[t], EA[t], a[t], e[t], b[t], p[t]}, t]


I then attempted this problem using ParametricNDSolve, and then subsequently fitting to the parameters using NonlinearModelFit and FindFit. With FindFit, I get the error:

"The function value {<<1>>} is not a list of real numbers with \
dimensions "


With NonlinearModelFit, I get

NDSolve Parametric Function is not a list of real \ numbers with dimensions


What I'm doing is here:

pSoln = ParametricNDSolve[
{EA'[t] ==
k1* e[t]* a[t] - k11* EA[t] + k22 *EAB[t] - k2* b[t] *EA[t] ,
EB'[t] == k4 e[t] b[t] - k44  EB[t] + k33 EAB[t] - k3 a[t] EB[t],
EAB'[t] ==
k2 EA[t] b[t] + k3 EB[t] a[t] - k22 EAB[t] - k5 EAB[t] +
k33 EAB[t],
p'[t] == k5 EAB[t],
e'[t] == k44 EB[t] + k11 EA[t] - k4 e[t] b[t] - k1 e[t] a[t] ,
a'[t] == k11 EA[t] - k1 e[t] a[t] + k33 EAB[t] - k3 EB[t] a[t],
b'[t] == k44 EB [t] - k4 e[t] b[t] + k22 EAB[t] - k2 EA[t] b[t] ,

EA[0] == 0,
EB[0] == 0,
EAB[0] == 0,
a[0] == 1,
e[0] == .1,
b[0] == 100,
p[0] == 0},
p, {t, 0, 1000}, {k1, k11, k2, k22, k3, k33, k4, k44, k5}]

fit = FindFit[data,
pSoln[k1, k11, k2, k22, k3, k33, k4, k44, k5][
t], {{k1, .15}, {k11, 1}, {k2}, k22, k3,
k33, {k4, .75}, {k44, 0}, {k5, .75}}, t];


For the NonlinearModelFit error, I just changed FindFit to NonlinearModelFit.

UPDATE: I think that ParametricNDSolve is the way to go based on the other stack posts I've viewed; however, I'm unable to plot the output of ParametricNDSolve/ParametricNDSolveValue. I had made the following function where I could tune the values of the constants:

    pFunction[t_?NumericQ, k1_?NumericQ, k11_?NumericQ, k2_?NumericQ,
k22_?NumericQ, k3_?NumericQ, k33_?NumericQ, k4_?NumericQ,
k44_?NumericQ, k5_?NumericQ]  :=
Evaluate[pSoln][k1, k11, k2, k22, k3, k33, k4, k44, k5][t];


With that, I then wanted to do NMinimize using the values from the manual sliders to find the actual fits:

{rmse, params} = NMinimize[{Sqrt[
Mean[(pFunction[data[[1]][[All, 1]], k1, k11, k2, k22, k3,
k33, k4, k44, k5] - data[[1]][[All, 2]])^2]],
0 < k1 < 0.2,
0 < k11 < 2,
0 < k2 < 0.05,
0 < k22 < 1,
0 < k3 < 0.2,
0 < k33 < 0.2,
0 < k4 < 2,
0 < k44 < 3,
0 < k5 < 1},
{k1, k11, k2, k22, k3, k33, k4, k44, k5}
]


Whether or not I used ?NumericQ or _Real for the inputs in pFunction, NMinimize crashes with

"The function value .... is \
not a number at {k1,k11,k2,k22,k3,k33,k4,k44,k5}"


PS. The question was edited.

• I should backtrack more here in the explanation to justify the approach that I am taking: I'm trying to fit a model to some data, where the model is a system of differential equations. So I wanted to get a symbolic solution so that I could do FindMinimum between my data and the model based on different sets of k values. All the stack posts I've read for the problem I post above with Solve suggest using some form of NDSolve. Using ParametricNDSolve, I'm able to get solutions very quickly, but I'm unsure how I can compare the output here with my actual data to backtrack the k constants. – wiscoYogi May 31 '19 at 21:40
• Is that eb[t] term supposed to be EB[t] or e[t] b[t]? – Carl Woll Jun 1 '19 at 0:05
• I've replaced this to EB[t]! No different in DSolve output though... – wiscoYogi Jun 1 '19 at 0:39
• Try using ParametricNDSolveValue instead, although with 9 parameters, I expect it will take a long time. – Carl Woll Jun 1 '19 at 1:45
• I've tried this (1) I don't know how to get the k values (2) it returns immediately – wiscoYogi Jun 1 '19 at 1:52

# Update

I was a little too clever when I extracted the dependent variables from the equations and made an error to the variable to fit. The 'p' variable to fit occurred in equation 4, but was extracted as the 7th variable. I have corrected the error so that the fit corresponds to the "p" variable.

As I commented in my answer, parameter fitting is frequent topic on MSE. I have been experimenting with StringTemplate as a possible way generalize the fitting to address this common pattern. Parameter fitting is a difficult subject and will depend on your data quality, your model, and your intial guesses. In your particular case, we are trying to fit 9 parameters to noisy data that could be well characterized by a straight line.

# Approach

• Use ParametricNDSolveValue to create the model.
• Use StringTemplates to handle lists of parameters and variables.
• Generate a Manipulate slider model to debug model and understand the effects of parameter changes.
• Transfer initial guesses from manipulate to perform a fit.

# Implementation

I commented the code so I hope it self explanatory. First assign the constants and prep the data.

(* Data for fitting Column 2 is "p" *)
(* data not reproduced here *)
(* Max Time *)
tmax = 1000;
(* Equations used to extract parameters and variables *)
test = {EA'[t] ==
k1*e[t]*a[t] - k11*EA[t] + k22*EAB[t] - k2*b[t]*EA[t],
EB'[t] == k4 e[t] b[t] - k44 EB[t] + k33 EAB[t] - k3 a[t] EB[t],
EAB'[t] ==
k2 EA[t] b[t] + k3 EB[t] a[t] - k22 EAB[t] - k5 EAB[t] +
k33 EAB[t], p'[t] == k5 EAB[t],
e'[t] == k44 EB[t] + k11 EA[t] - k4 e[t] b[t] - k1 e[t] a[t],
a'[t] == k11 EA[t] - k1 e[t] a[t] + k33 EAB[t] - k3 EB[t] a[t],
b'[t] == k44 EB[t] - k4 e[t] b[t] + k22 EAB[t] - k2 EA[t] b[t]};


Here is the code to solve the Diffeqs and create the Manipulate slider to improve initial guesses.

(* Generate System of Differential Equations *)
e1 = EA'[t] ==
k1*e[t]*a[t] - k11*EA[t] + k22*EAB[t] - k2*b[t]*EA[t]; e2 =
EB'[t] == k4 e[t] b[t] - k44 EB[t] + k33 EAB[t] - k3 a[t] EB[t]; e3 =
EAB'[t] ==
k2 EA[t] b[t] + k3 EB[t] a[t] - k22 EAB[t] - k5 EAB[t] +
k33 EAB[t]; e4 = p'[t] == k5 EAB[t];
e5 = e'[t] ==
k44 EB[t] + k11 EA[t] - k4 e[t] b[t] - k1 e[t] a[t]; e6 =
a'[t] == k11 EA[t] - k1 e[t] a[t] + k33 EAB[t] - k3 EB[t] a[t]; e7 =
b'[t] == k44 EB[t] - k4 e[t] b[t] + k22 EAB[t] - k2 EA[t] b[t];
ics = {EA[0] == 0, EB[0] == 0, EAB[0] == 0, a[0] == 1, e[0] == 1/10,
b[0] == 100, p[0] == 0};
eqns = {e1, e2, e3, e4, e5, e6, e7}~Join~ics;
(* Extract Variables from test *)
vbles = Cases[test, f_'[t_] -> f, All] // Union;
(*Parameters with target and desired ranges*)
mat = ConstantArray[0, {9, 4}];
(* Extract parameters from test *)
mat[[All, 1]] = Cases[test, _Symbol, Infinity] // Union // Most;
mat = {
{k1, 0.15, 0.01, 0.4},
{k11, 1, 0.01, 2},
{k2, 0.05, 0, 0.1},
{k22, 0.05, 0, 0.1},
{k3, 0.05, 0, 0.1},
{k33, 0.05, 0, 0.5},
{k4, 0.75, 0.01, 1},
{k44, 0.5, 0, 1},
{k5, 0.75, 0, 2}
};
(* Reduced parameter matrix if we want to fix variables *)
rmat = mat;
(* Build Manipulate sliders *)
sfun =  StringRiffle[(StringTemplate[
"{{1,2},3,4,Appearance\[Rule]\"Labeled\"}"] @@ #) & \
/@ #, ","] &;
sliders = sfun[rmat];
(* Extract Parameters from mat *)
parms = mat[[All, 1]];
(* Create String Representations of parms *)
sparms = StringRiffle[ToString[#] & /@ parms, ","];
(* Create patterns and string reps of parameters *)
pats = Pattern @@@ (#*_ & /@ parms);
spats = StringRiffle[ToString[#] & /@ pats, ","];
(* List Plot of the data *)
lp = ListPlot[data, PlotLegends -> {"data-p"}];
(* ParametricNDSolveValue *)
pfun = ParametricNDSolveValue[eqns, vbles, {t, 0, tmax}, parms];
(*Create an appropriate model function to fit*)
modelstring = "{#[[1]],#[[2]],#[[3]],#[[4]],#[[5]],#[[6]],#[[7]]}&";
(* Create some PlotLegends *)
pl = ",PlotLegends\[Rule]" <> StringRiffle[ToString[#], ","] &@vbles;
(* Build the model expression *)
ToExpression[
StringTemplate[
"model[pats][t_]:=ms@Through[pfun[params][t],List]\
/;And@@NumericQ/@{params};"][<|"pats" -> spats, "params" -> sparms,
"ms" -> modelstring|>]]
(* Create slider model to help with initial guess *)
globalstring =
StringTemplate["global={params};"][<|"params" -> sparms|>];
mantemp =
"Manipulate[g\[IndentingNewLine]Show[lp,Plot[Evaluate@(model[\
params][t]),{t,0,tmax},PlotRange\[Rule]{0.,.8}pl],ImageSize->Large]\
,sliders]";
ToExpression@
StringTemplate[mantemp][<|"sliders" -> sliders, "params" -> sparms,
"g" -> globalstring, "pl" -> pl|>]
(*Display global variable of current slider parameters *)
Dynamic@global


We can adjust the sliders to get a better initial guess for the parameters by trying to manually the sliders from the initial values. I intentionally did not manually fit as well as I could have to test how well FindFit change parameters

Here is the code to do the fit and create a slider model initialized at the starting parameters. The "Gradient" method seemed to work reasonably well, but the methods will need to be experimented with to achieve the best results.

(*Initial parameter guesses from global*)
fit = Quiet@
FindFit[data,
model[k1, k11, k2, k22, k3, k33, k4, k44, k5][t][[7]], initguess,
(*Create new Manipulate Based on Fit*)
mat2 = mat;
mat2[[All, 2]] = fit[[All, 2]];
ToExpression@
StringTemplate[mantemp][<|"sliders" -> sfun[mat2],
"params" -> sparms, "pl" -> pl|>]


I perturbed the initial guesses and re-ran the fit. The fit looks similar, but the fitted parameters are clearly different indicating the non-uniqueness depending on initial guesses. Good initial and long term asymptotic data would aid in the parameter estimation.

• Can you clarify how model can be called while fitting/validating the starting model with the constants determined from the initial sliders in the part above (above the second image)? Syntactically I understand that model leverages the output of ParametricNDSolveValue as well as the guesses of the constants, but I ask because model is defined in strings in the first part, and so when typing model and pressing shift + return, nothing happens. Thank you so much for this thorough response. – wiscoYogi Jun 1 '19 at 18:10
• ^can't edit this... calling model[[4]] results in the error message "part 4 of model does not exist"... I can achieve this with model[[1]] – wiscoYogi Jun 1 '19 at 18:43
• I will have a look. Sometimes the text can get garbled a little bit when copied to MSE. You can see the current definition of model by ?model. To convert the model string to an expression you need to issue ToExpression[ StringTemplate[ "model[pats][t_]:=ms@Through[pfun[params][t],List]\ /;And@@NumericQ/@{params};"][<|"pats" -> spats, "params" -> sparms, "ms" -> modelstring|>]]. – Tim Laska Jun 1 '19 at 20:18
• I also didn't understand why FindFit is the best way to make agreement with the model. I tried using NMinimize, will put into the post above. – wiscoYogi Jun 1 '19 at 20:19
• FindFit is supposed to be a simple wrapper to other more complex optimizers. You can try Method->NMinimize with FindFit and it will invoke the default NMinimize methods. It runs forever on my machine. I also think that we should be using part 7 "variable p" versus 4 "variable EA". Can you confirm that your data is p vs t? – Tim Laska Jun 1 '19 at 21:26

You can get a series expansion for small t with AsymptoticDSolveValue; maybe that's good enough if you go to high-enough order in t:

AsymptoticDSolveValue[{
EA'[t] == k1 e[t] a[t] + k22 EAB[t] - k2 b[t] EA[t] - k11 EA[t],
EB'[t] == k4 e[t] b[t] - k44 EB[t] + k33 EAB[t] - k3 a[t] eb[t],
EAB'[t] == k2 EA[t] b[t] + k3 EB[t] a[t] - k22 EAB[t] - k5 EAB[t] + k33 EAB[t],
p'[t] == k5 EAB[t],
e'[t] == k44 EB[t] + k11 EA[t] - k4 e[t] b[t] - k1 e[t] a[t],
a'[t] == k11 EA[t] - k1 e[t] a[t] + k33 EAB[t] - k3 EB[t] a[t],
b'[t] == k44 EB[t] - k4 e[t] b[t] + k22 EAB[t] - k2 EA[t] b[t],
EA[0] == 0, EB[0] == 0, EAB[0] == 0, a[0] == 1, e[0] == 1/10, b[0] == 100, p[0] == 0},
{EAB[t], EB[t], EA[t], a[t], e[t], b[t], p[t]},
{t, 0, 3}]


(output is third-order series expansion in t)

• I need to go up to t=1000, so this isn't going to work... I'm trying to use this approach with the Parametric NDSolve, but I keep getting "the function value is not a list of real numbers with dimensions" community.wolfram.com/groups/-/m/t/126143 – wiscoYogi May 31 '19 at 22:06
• Yes I agree that with your specific parameters this asymptotic expansion does not converge. In other contexts it's a handy tool though. – Roman Jun 2 '19 at 19:57