1
$\begingroup$

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?

$\endgroup$
6
  • $\begingroup$ your Manipulate is incomplete. $\endgroup$
    – george2079
    Nov 6, 2017 at 15:20
  • $\begingroup$ 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. $\endgroup$
    – m_goldberg
    Nov 6, 2017 at 15:39
  • $\begingroup$ 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. $\endgroup$
    – OKCarl
    Nov 6, 2017 at 17:17
  • $\begingroup$ Without you providing the points (e.g. by posting them on Pastebin), no one can help you with your problem. $\endgroup$ Nov 6, 2017 at 21:09
  • $\begingroup$ 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 $\endgroup$
    – OKCarl
    Nov 6, 2017 at 22:41

1 Answer 1

1
$\begingroup$

You can use MCMC package which can be downloaded here https://github.com/stan-dev/MathematicaStan

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]]

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.