# Problem using NonlinearModelFit to 2 coupled differential equations

I am trying to find the parameters of a set of differential equations given experimental data. Before I tackled my set of 12 equations, I tried to solve a simpler and similar 2 equation system but have not been able to manage.

First I created a system and solved it, summed the outputs{x,y} and then added noise.

Clear[f, g, x, y, model]
sol = NDSolve[{y''[t] == 5*x[t] - y[t], x''[t] == 2*y[t] - 10*x[t],
y == 3, y' == 0, x == 3, x' == 5}, {x, y}, {t, 0,
5}];

tvals = Range[0, 5, 0.1];

data = Transpose[{tvals,
Flatten[Evaluate[x[tvals] /. sol] + Evaluate[y[tvals] /. sol]] +
RandomReal[{-0.5, 0.5}, Length[tvals]]}];

Show[Plot[Evaluate[y[t] /. sol] + Evaluate[x[t] /. sol], {t, 0, 5},
PlotStyle -> Green], ListPlot[data, PlotRange -> All]]


Next, I try to fit a curve and find the parameters of the original differential equations {5,-1,2,-10}. Here is my attempt:

model[a_?NumberQ, b_?NumberQ, c_?NumberQ,
d_?NumberQ] := (model[a, b, c, d] =
First[sol2 =
NDSolve[{f''[t] == a*g[t] - b*f[t], g''[t] == c*f[t] - d*g[t],
f == 3, f' == 0, g == 3, g' == 5},
g + f, {t, 0, 5}]])

nlm = NonlinearModelFit[data,
Evaluate[model[a, b, c, d][t] /. sol2], {a, b, c, d}, t];

Show[Plot[nlm[t], {t, 0, 5}, PlotRange -> All], ListPlot[data]]
nlm["ParameterTable"]


All I get is error messages an I am unsure of what to do.

• Your model does not depend on a,b,c,d... – Marius Ladegård Meyer Aug 7 '18 at 16:26
• Also, as long as you have the _?NumberQ pattern, the Evaluate in NonlinearModelFit won't do anything, because a,b,c,d won't be numbers when Evaluate is evoked. – Marius Ladegård Meyer Aug 7 '18 at 16:27

In this answer I assume that the variables {a, b, c, d} are given in the differential equations model. (Not as initial values as in this answer.)

The idea is to make set of functions to fit with different parameters and pick the ones that give the best fit.

Redefine the model:

Clear[model]
model[a_?NumberQ, b_?NumberQ, c_?NumberQ, d_?NumberQ] :=
First[sol2 =
NDSolve[{f''[t] == b*g[t] - a *f[t], g''[t] == c*f[t] - d*g[t],
f == 3, f' == 0, g == 3, g' == 5}, {f, g}, {t, 0,
5}]];

model[1, 2, 3, 2]


Construct a (reasonable) set of functions:

AbsoluteTiming[
bFuncs =
Association[
Flatten[Table[{a, b, c, d} -> model[a, b, c, d], {a, 0, 6, 1}, {b, 0, 10, 1}, {c, 0, 6, 1}, {d, 8, 12, 1}]]];
]
(* {3.31016, Null} *(

Length[bFuncs]

(* 2695 *)


Compute the fit for all functions in the set:

AbsoluteTiming[
fitErrRes =
Table[(
fit = Fit[data, Through[Values[bFuncs[[i]]][t]], t];
errs = Map[(fit /. t -> #[]) - #[] &, data];
Keys[bFuncs][[i]] -> Norm[errs, \[Infinity]]), {i, 1,
Length[bFuncs]}];
]

(* {5.37995, Null} *)


Pick the bet fit:

Through[{Min, Max, Mean, StandardDeviation}[
Values[fitErrRes]]]

(* {0.538426, 23.9603, 15.4338, 4.54161} *)

sInd =
Flatten[Position[Values[fitErrRes], Min[Values[fitErrRes]]]][]

(* 573 *)

Keys[bFuncs][[sInd]]

(* {1, 5, 2, 10} *)


Plot

Show[{ListPlot[data],
Plot[Evaluate[Total[Through[Values[bFuncs[[sInd]]][t]]]], {t,
Min[data[[All, 1]]], Max[data[[All, 1]]]}]}] You can extend and/or narrow the search by changing the set search functions bFuncs.

• Thanks all. A mistake on my behalf. The {a,b,c,d} were meant to be the coefficients in the DEs but I was playing around with the code and forgot to change the numbers back to letters before I copied and pasted the code. – Dave Bassito Aug 8 '18 at 12:26
• @VincenzoPratley It is a good question. If your code was more streamlined I assume it would have got more attention. Please make the corresponding clarifications... – Anton Antonov Aug 8 '18 at 12:58
• edits made. I have only been using Mathematica for about 2 weeks and I am also a relative beginner at programming. – Dave Bassito Aug 8 '18 at 14:50
• @VincenzoPratley Ok, thank you for doing that. – Anton Antonov Aug 8 '18 at 16:36
Clear[f, g, x, y, model]
sol = NDSolve[{y''[t] == 5*x[t] - y[t], x''[t] == 2*y[t] - 10*x[t],
y == 3, y' == 0, x == 3, x' == 5}, {x, y}, {t, 0, 5}];

tvals = Range[0, 5, 0.1];

data = Transpose[{tvals,
Flatten[Evaluate[x[tvals] /. sol] + Evaluate[y[tvals] /. sol]] +
RandomReal[{-0.5, 0.5}, Length[tvals]]}];

Show[Plot[Evaluate[y[t] /. sol] + Evaluate[x[t] /. sol], {t, 0, 5},
PlotStyle -> Green], ListPlot[data, PlotRange -> All]]
model[a_?NumberQ, b_?NumberQ, c_?NumberQ, d_?NumberQ] :=
Module[{y, t},
First[y /.
NDSolve[{f''[t] == 5*g[t] - f[t], g''[t] == 2*f[t] - 10*g[t],
f == a, f' == b, g == c, g' == d,
y[t] == f[t] + g[t], y == a + c}, {g, f, y}, {t, 0, 5}]]] //
Quiet

nlm = NonlinearModelFit[data, model[a, b, c, d][t], {a, b, c, d},
t]; // Quiet

Show[Plot[nlm[t], {t, 0, 5}, PlotRange -> All], ListPlot[data]]
nlm["ParameterTable"] • Good (+1). Are was thinking that {a,b,c,d} are coefficients in the differential equations, not in the initial values. – Anton Antonov Aug 7 '18 at 18:20