# NonlinearModelFit given the solution of a differential equation as its model

I have a complicated differential equation working quite well in a Manipulate expression. Here is my code:

ClearAll[t, E6a, x0, v0, a1, f1, m1, c1, k0, k1, e1]

E6a[t_, x0_, v0_, a1_, f1_, m1_, c1_, k0_, k1_, e1_] :=
xx[t] /.
NDSolve[
{m1 xx''[t] + c1 xx'[t] + (k1 + (k0 - k1) Exp[-t e1]) xx[t] ==
a1 Sin[f1 2 π ( t/125)],
xx[0] == x0, xx'[0] == v0},
xx, {t, 0, 250}]


Unfortunately, although I have good start values for the variables, I am not able to get a reasonable solution for the parameters from NonlinearModelFit. Here is the code I used for that:

Import["https://pastebin.com/raw/k4ume1yR", "Package"];

varElast =
NonlinearModelFit[
ptsM50thIn50th,
E6a[x0, v0, a1, f1, m1, c1, k0, k1, e1][t],
{{x0, -20}, {v0, -5700}, {a1, 10.0}, {f1, 0.725}, {m1, 0.01},
{c1, 9.0}, {k0, 3.25}, {k1, 1.5}, {e1, 3.3}},
t]


Could someone please tell me, what am I doing wrong? Where is my mistake?

• your Manipulate is incomplete. Commented Nov 6, 2017 at 15:20
• 1) I have removed your incomplete Manipulate code because I could not see its relevance to your question. 2) I don't see how to help you without having a definition of ptsM50thIn50th. Commented Nov 6, 2017 at 15:39
• Thank you for taking the time to read my post. In answer to your question: "ptsM50thIn50th" is a set of {X, Y} data points. As for example, my first set of data points = {9.5, -18.93} a later set of points: {105, 0.666}. Thanks again. Commented Nov 6, 2017 at 17:17
• Without you providing the points (e.g. by posting them on Pastebin), no one can help you with your problem. Commented Nov 6, 2017 at 21:09
• Dear JM, Thanks for the info. I posted my points on Pastebin under the title: OKCarl's Data Point Set. Here is (I think) the direct link: pastebin.com/k4ume1yR Commented Nov 6, 2017 at 22:41

put the .nb and .m file in the same folder. It can be improved by choosing different initial points. Good luck.

Get[FileNameJoin[{NotebookDirectory[], "mcmc.m"}]]

data = {{9.5, -18.93}, {11, -13.14}, {13, -8.830}, {15, -5.940}, {17, \
-3.758}, {19, -2.035}, {21, -0.6384}, {23, 0.5206}, {25, 1.5027}, {27,
2.3182}, {29, 3.0175}, {31, 3.601}, {33.5, 3.405}, {37,
3.989}, {41, 4.381}, {45, 4.537}, {49, 4.572}, {53, 4.507}, {57,
4.346}, {61, 4.112}, {65, 3.828}, {69, 3.527}, {73, 3.213}, {77,
2.888}, {81, 2.566}, {85, 2.259}, {89, 1.951}, {93, 1.644}, {97,
1.330}, {101, 1.004}, {105, 0.666}, {109,
0.310}, {113, -0.083}, {117, -0.490}, {121, -0.930}, {125, \
-1.380}, {129, -1.850}, {133, -2.297}, {137, -2.710}, {141, -3.049}, \
{145, -3.282}, {149, -3.369}, {153, -3.259}, {157, -2.921}, {161, \
-2.345}, {165, -1.551}, {169, -0.599}, {173, 0.422}, {177,
1.390}, {181, 2.227}, {185, 2.841}, {189, 3.202}, {193,
3.319}, {197, 3.192}, {201, 2.881}, {205, 2.414}, {209,
1.821}, {213, 1.140}, {217,
0.402}, {221, -0.376}, {225, -1.173}, {229, -1.972}, {233, \
-2.772}, {237, -3.566}, {241, -4.343}, {245, -5.106}, {249, -5.845}, \
{249.5, -5.933}};

model[a1_?NumberQ, f1_?NumberQ, m1_?NumberQ, c1_?NumberQ, k0_?NumberQ,
k1_?NumberQ,
e1_?NumberQ] := (model[a1, f1, m1, c1, k0, k1, e1] =
Module[{x, t},
NDSolveValue[{m1 x''[t] +
c1 x'[t] + (k1 + (k0 - k1) Exp[-t e1]) x[t] ==
a1 Sin[f1 2 \[Pi]/125 t], x[0] == -20, x'[0] == -5700},
x, {t, 0, 250}]]);

error = ConstantArray[0.5, Length@data];

spr = 0.01; mcmc =
MCMCModelFit[data, error,
model[a1, f1, m1, c1, k0, k1, e1][
t], {{a1, 10, spr, Reals}, {f1, 0.7, spr, Reals}, {m1, 0.01, spr,
Reals}, {c1, 9, spr, Reals}, {k0, 3, spr, Reals}, {k1, 1.5, spr,
Reals}, {e1, 3, spr, Reals}}, {t}, 10000];
mcmc["BestFitParameters"]

{a1 -> 7.14169, f1 -> 0.830792, m1 -> 0.223974, c1 -> 11.3987,
k0 -> 2.69359, k1 -> 2.23043, e1 -> 3.24765}

Show[Plot[
model[a1, f1, m1, c1, k0, k1, e1][t] /.
mcmc["BestFitParameters"], {t, 0, 250}, PlotTheme -> "Detailed",
ImageSize -> 400, PlotLegends -> {"model"},
PlotRange -> {MinMax@data[[All, 1]],
MinMax@data[[All, 2]] + {-1, 1}}],
ListPlot[data, PlotStyle -> Red]]