4
$\begingroup$

I have this differential equation: $$m\ddot x=-kx^\frac{3}{2}-c\dot x-mg$$ where I want to fit for $k$, $c$. ($g$ is 9.81 and $m$ is 0.3).

This is a model for collision, hence in the data that we have collected in our experiment, all we know is that x'[0]==-3, where -3 is the impact velocity before the collision, and x'[T]==2 where 2 is the rebounding velocity after the collision and Tis the time of contact, which we cannot measure experimentally as it is very short, but we know that it is shorter than $10^{-3}s$.

m = 1;
k = 1;
c = 1;
g = 9.81;
sol = NDSolve[ 
  {m x''[t] == -k x[t]^(3/2) - c x'[t] - m g, x'[0] == -3, x[0] == 0.024965, 
   x'[0.00001] == 2},
  x[t], {t, 0, 1}]

Here is the data.

Data for x against t:

{{0.,23.6724},{0.0333333,23.4316},{0.0666667,23.2125},
 {0.1,22.9737},{0.133333,22.7191},{0.166667,22.4796},
 {0.2,22.2635},{0.233333,22.0175},{0.266667,21.7774},
 {0.3,21.5224},{0.333333,21.3139},{0.366667,21.064},
 {0.4,20.8183},{0.433333,20.5699},{0.466667,20.3129},
 {0.5,20.0644},{0.533333,19.8333},{0.566656,19.5862},
 {0.599989,19.3391},{0.633322,19.094},{0.666656,18.8495},
 {0.699989,18.5973},{0.733322,18.3451},{0.766656,18.09},
 {0.799989,17.8299},{0.833322,17.581},{0.866656,17.3204},
 {0.899989,17.0659},{0.933322,16.817},{0.966656,16.5627},
 {0.999989,16.3046},{1.03332,16.0535},{1.06666,15.7956},
 {1.09999,15.5383},{1.13332,15.2806},{1.16666,15.0236},
 {1.19999,14.7635},{1.23332,14.5015},{1.26666,14.2514},
 {1.29999,13.9673},{1.33332,13.6998},{1.36666,13.4402},
 {1.39999,13.1574},{1.43332,12.8848},{1.46666,12.6188},
 {1.49999,12.3376},{1.53332,12.0596},{1.56666,11.7867},
 {1.59999,11.5302},{1.63332,11.2418},{1.66664,10.9721},
 {1.69998,10.7005},{1.73331,10.399},{1.76664,10.1111},
 {1.79998,9.83385},{1.83331,9.56173},{1.86664,9.25114},
 {1.89998,8.98928},{1.93331,8.70041},{1.96664,8.41822},
 {1.99998,8.13319},{2.03331,7.84509},{2.06664,7.53343},
 {2.09998,7.25237},{2.13331,6.95413},{2.16664,6.63875},
 {2.19998,6.34642},{2.23331,6.06828},{2.26664,5.77579},
 {2.29998,5.4747},{2.33331,5.15976},{2.36664,4.84916},
 {2.39998,4.5256},{2.43331,4.22336},{2.46664,3.9177},
 {2.49998,3.58284},{2.53331,3.2908},{2.56664,2.97411},
 {2.59998,2.6861},{2.63331,2.4965},{2.66664,2.73492},
 {2.69998,2.99366},{2.73331,3.29602},{2.76663,3.58096},
 {2.79997,3.83507},{2.8333,4.1179},{2.86663,4.39381},
 {2.89997,4.66047},{2.9333,4.95059},{2.96663,5.23038},
 {2.99997,5.48554},{3.0333,5.77507},{3.06663,6.03556},
 {3.09997,6.30288},{3.1333,6.56806},{3.16663,6.82612},
 {3.19997,7.11681},{3.2333,7.37396},{3.26663,7.63213},
 {3.29997,7.89755},{3.3333,8.15167},{3.36663,8.4428},
 {3.39997,8.6969},{3.4333,8.95516},{3.46663,9.22325},
 {3.49997,9.47407},{3.5333,9.73972},{3.56663,9.98549},
 {3.59997,10.2457},{3.6333,10.4917},{3.66663,10.7494},
 {3.69997,10.9985},{3.7333,11.2493},{3.76663,11.5069},
 {3.79997,11.7599},{3.8333,12.0148},{3.86663,12.2645},
 {3.89996,12.5198},{3.93329,12.7714},{3.96662,13.0222},
 {3.99996,13.2753},{4.03329,13.4973},{4.06662,13.7457},
 {4.09996,13.9856},{4.13329,14.2364},{4.16662,14.4828},
 {4.19996,14.7348},{4.23329,14.9753},{4.26662,15.211},
 {4.29996,15.4466},{4.33329,15.6922},{4.36662,15.9198},
 {4.39996,16.1627},{4.43329,16.4001},{4.46662,16.6353},
 {4.49996,16.8629},{4.53329,17.1011},{4.56662,17.3418},
 {4.59996,17.5674},{4.63329,17.81},{4.66662,18.0313},
 {4.69996,18.2533},{4.73329,18.4823},{4.76662,18.7227},
 {4.79996,18.9488},{4.83329,19.1835},{4.86662,19.4019},
 {4.89996,19.6282},{4.93329,19.86},{4.96662,20.084},
 {4.99994,20.3083},{5.03328,20.5353},{5.06661,20.7602},
 {5.09994,20.9745},{5.13328,21.1844},{5.16661,21.4296},
 {5.19994,21.6461},{5.23328,21.8579},{5.26661,22.0885},
 {5.29994,22.3081},{5.33328,22.5211}}

Take note that x is in cm.

Most of the data is useless because they are just data for the dropping and bouncing part, not actually the collision.

In the code, I only did NDSolve and substitute in random values for $k$, $c$, and also substitute some of the initial conditions like x[0]==0.024965 , x'[0]==-3 and x[T]==2.

With these, is it possible for us to fit the constants?

Thank you.

$\endgroup$
14
  • $\begingroup$ Please post the data you have collected in your experiment. $\endgroup$ – Anton Antonov Dec 10 '20 at 1:46
  • $\begingroup$ Okay sure I will edit the post. But i dont think the data I have collected would be helpful. Because the part that actually matters only take up less than one frame $\endgroup$ – bob the legend Dec 10 '20 at 1:47
  • $\begingroup$ Hi sir does it work? $\endgroup$ – bob the legend Dec 10 '20 at 3:58
  • $\begingroup$ It can be made to work (to a point), but multiple adjustments have to be made. I will post an answer tomorrow. $\endgroup$ – Anton Antonov Dec 10 '20 at 5:14
  • $\begingroup$ Thank you so much. $\endgroup$ – bob the legend Dec 10 '20 at 5:31
2
$\begingroup$

Actually we can use data to optimise parameters as follows

data = {{0., 23.6724}, {0.0333333, 23.4316}, {0.0666667, 23.2125}, {0.1, 22.9737}, {0.133333, 22.7191}, {0.166667, 22.4796}, {0.2, 22.2635}, {0.233333, 22.0175}, {0.266667, 21.7774}, {0.3, 21.5224}, {0.333333, 21.3139}, {0.366667, 21.064}, {0.4, 20.8183}, {0.433333, 20.5699}, {0.466667, 20.3129}, {0.5, 20.0644}, {0.533333, 19.8333}, {0.566656, 19.5862}, {0.599989, 19.3391}, {0.633322, 19.094}, {0.666656, 18.8495}, {0.699989, 18.5973}, {0.733322, 18.3451}, {0.766656, 18.09}, {0.799989, 17.8299}, {0.833322, 17.581}, {0.866656, 17.3204}, {0.899989, 17.0659}, {0.933322, 16.817}, {0.966656, 16.5627}, {0.999989, 16.3046}, {1.03332, 16.0535}, {1.06666, 15.7956}, {1.09999, 15.5383}, {1.13332, 15.2806}, {1.16666, 15.0236}, {1.19999, 14.7635}, {1.23332, 14.5015}, {1.26666, 14.2514}, {1.29999, 13.9673}, {1.33332, 13.6998}, {1.36666, 13.4402}, {1.39999, 13.1574}, {1.43332, 12.8848}, {1.46666, 12.6188}, {1.49999, 12.3376}, {1.53332, 12.0596}, {1.56666, 11.7867}, {1.59999, 11.5302}, {1.63332, 11.2418}, {1.66664, 10.9721}, {1.69998, 10.7005}, {1.73331, 10.399}, {1.76664, 10.1111}, {1.79998, 9.83385}, {1.83331, 9.56173}, {1.86664, 9.25114}, {1.89998, 8.98928}, {1.93331, 8.70041}, {1.96664, 8.41822}, {1.99998, 8.13319}, {2.03331, 7.84509}, {2.06664, 7.53343}, {2.09998, 7.25237}, {2.13331, 6.95413}, {2.16664, 6.63875}, {2.19998, 6.34642}, {2.23331, 6.06828}, {2.26664, 5.77579}, {2.29998, 5.4747}, {2.33331, 5.15976}, {2.36664, 4.84916}, {2.39998, 4.5256}, {2.43331, 4.22336}, {2.46664, 3.9177}, {2.49998, 3.58284}, {2.53331, 3.2908}, {2.56664, 2.97411}, {2.59998, 2.6861}, {2.63331, 2.4965}, {2.66664, 2.73492}, {2.69998, 2.99366}, {2.73331, 3.29602}, {2.76663, 3.58096}, {2.79997, 3.83507}, {2.8333, 4.1179}, {2.86663, 4.39381}, {2.89997, 4.66047}, {2.9333, 4.95059}, {2.96663, 5.23038}, {2.99997, 5.48554}, {3.0333, 5.77507}, {3.06663, 6.03556}, {3.09997, 6.30288}, {3.1333, 6.56806}, {3.16663, 6.82612}, {3.19997, 7.11681}, {3.2333, 7.37396}, {3.26663, 7.63213}, {3.29997, 7.89755}, {3.3333, 8.15167}, {3.36663, 8.4428}, {3.39997, 8.6969}, {3.4333, 8.95516}, {3.46663, 9.22325}, {3.49997, 9.47407}, {3.5333, 9.73972}, {3.56663, 9.98549}, {3.59997, 10.2457}, {3.6333, 10.4917}, {3.66663, 10.7494}, {3.69997, 10.9985}, {3.7333, 11.2493}, {3.76663, 11.5069}, {3.79997, 11.7599}, {3.8333, 12.0148}, {3.86663, 12.2645}, {3.89996, 12.5198}, {3.93329, 12.7714}, {3.96662, 13.0222}, {3.99996, 13.2753}, {4.03329, 13.4973}, {4.06662, 13.7457}, {4.09996, 13.9856}, {4.13329, 14.2364}, {4.16662, 14.4828}, {4.19996, 14.7348}, {4.23329, 14.9753}, {4.26662, 15.211}, {4.29996, 15.4466}, {4.33329, 15.6922}, {4.36662, 15.9198}, {4.39996, 16.1627}, {4.43329, 16.4001}, {4.46662, 16.6353}, {4.49996, 16.8629}, {4.53329, 17.1011}, {4.56662, 17.3418}, {4.59996, 17.5674}, {4.63329, 17.81}, {4.66662, 18.0313}, {4.69996, 18.2533}, {4.73329, 18.4823}, {4.76662, 18.7227}, {4.79996, 18.9488}, {4.83329, 19.1835}, {4.86662, 19.4019}, {4.89996, 19.6282}, {4.93329, 19.86}, {4.96662, 20.084}, {4.99994, 20.3083}, {5.03328, 20.5353}, {5.06661, 20.7602}, {5.09994, 20.9745}, {5.13328, 21.1844}, {5.16661, 21.4296}, {5.19994, 21.6461}, {5.23328, 21.8579}, {5.26661, 22.0885}, {5.29994, 22.3081}, {5.33328, 22.5211}};

Now we can use interpolation function f = Interpolation[data, InterpolationOrder -> 4] to find out dependence of acceleration on x and x' as

{ParametricPlot[{f[t], f''[t]}, {t, 2.55, 2.7}, PlotRange -> All, 
  AspectRatio -> 1/2, AxesLabel -> {"x", "x''"}], 
 ParametricPlot[{f'[t], f''[t]}, {t, 2.3, 2.8}, PlotRange -> All, 
  AspectRatio -> 1/2, AxesLabel -> {"x'", "x''"}]} 

Figure 1

It looks like typical elastic-plastic deformation, and therefore the Hertz model is not applicable at all. Now we can propose force before and after collision in a form $$F/m=-k_1 x+k_2 x^2 + k_3 \dot {x}+k_4 \dot {x}^2-g $$ Finally, using f[t] we can optimize model in several points, for example,

g=981.; param = Table[{t, 
   NMinimize[{(f''[t] + g - k1 f[t] + k2 f[t]^2 + k3 f'[t] + 
        k4 f'[t]^2)^2, k1 > 0 && k2 > 0 && k3 > 0 && k4 > 0}, {k1, k2,
      k3, k4}]}, {t, 2.51, 2.7, .01}]

From this table we see that parameters of the model drastically change after collision at t=2.63

{ListLinePlot[
  Table[{param[[i, 1]], k1 /. param[[i, 2, 2]]}, {i, Length[param]}], 
  AxesLabel -> {"t", "k1"}], 
 ListLinePlot[
  Table[{param[[i, 1]], k2 /. param[[i, 2, 2]]}, {i, Length[param]}], 
  AxesLabel -> {"t", "k2"}], 
 ListLinePlot[
  Table[{param[[i, 1]], k3 /. param[[i, 2, 2]]}, {i, Length[param]}], 
  AxesLabel -> {"t", "k3"}], 
 ListLinePlot[
  Table[{param[[i, 1]], k4 /. param[[i, 2, 2]]}, {i, Length[param]}], 
  AxesLabel -> {"t", "k4"}, PlotRange -> All]}

Figure 2

$\endgroup$
12
  • $\begingroup$ Very nice/interesting! $\endgroup$ – Anton Antonov Dec 10 '20 at 16:54
  • $\begingroup$ Thanks this is really helpful! :) $\endgroup$ – bob the legend Dec 11 '20 at 0:03
  • $\begingroup$ I have a few question: 1. Why is hertz contact model not applicable here? 2. I run the same code in mathematica, (just copy pasted the code in your answer), but the t-axis of my graphs of k against t is extremely small, and graphs of k1,k2,and k4 are all empty. May I know why is that so? $\endgroup$ – bob the legend Dec 11 '20 at 2:31
  • $\begingroup$ @bobthelegend 1) the Hertz model describes elastic collision, and in this case we have elastic-plastic deformation (see functions x''[x], x''[x']; 2) what version do you run? Let try to evaluate param // TableForm and check first line in the table. Did you see k1 -> 241.097 k2 -> 185.712 k3 -> 1048.46 k4 -> 89.1876 ? $\endgroup$ – Alex Trounev Dec 11 '20 at 9:53
  • $\begingroup$ No... It outputs: {8.27181*10^-25, {k1 -> 447.611, k2 -> 317.74, k3 -> 1201.73, k4 -> 104.398}}. All I did was copy paste the code included in your answer, and added a g=9.81; at the front. what is wrong here? (´・ω・`)? $\endgroup$ – bob the legend Dec 11 '20 at 12:26
3
$\begingroup$

I know I'm a little bit late, but I want to show how to solve the physical problem straighforward, based on the measurement tx (in units s,m !)

tx = Map[{#[[1]], #[[2]]/100} &,
{{0., 23.6724}, {0.0333333,23.4316}, {0.0666667, 23.2125}, {0.1, 22.9737}, {0.133333, 22.7191}, {0.166667, 22.4796}, {0.2, 22.2635}, {0.233333,22.0175}, {0.266667, 21.7774}, {0.3, 21.5224}, {0.333333,21.3139}, {0.366667, 21.064}, {0.4, 20.8183}, {0.433333,20.5699}, {0.466667, 20.3129}, {0.5, 20.0644}, {0.533333,19.8333}, {0.566656, 19.5862}, {0.599989, 19.3391}, {0.633322,19.094}, {0.666656, 18.8495}, {0.699989, 18.5973}, {0.733322,18.3451}, {0.766656, 18.09}, {0.799989, 17.8299}, {0.833322,17.581}, {0.866656, 17.3204}, {0.899989, 17.0659}, {0.933322,16.817}, {0.966656, 16.5627}, {0.999989, 16.3046}, {1.03332,16.0535}, {1.06666, 15.7956}, {1.09999, 15.5383}, {1.13332,15.2806}, {1.16666, 15.0236}, {1.19999, 14.7635}, {1.23332,14.5015}, {1.26666, 14.2514}, {1.29999, 13.9673}, {1.33332,13.6998}, {1.36666, 13.4402}, {1.39999, 13.1574}, {1.43332,12.8848}, {1.46666, 12.6188}, {1.49999, 12.3376}, {1.53332,12.0596}, {1.56666, 11.7867}, {1.59999, 11.5302}, {1.63332,11.2418}, {1.66664, 10.9721}, {1.69998, 10.7005}, {1.73331,10.399}, {1.76664, 10.1111}, {1.79998, 9.83385}, {1.83331,9.56173}, {1.86664, 9.25114}, {1.89998, 8.98928}, {1.93331,8.70041}, {1.96664, 8.41822}, {1.99998, 8.13319}, {2.03331,7.84509}, {2.06664, 7.53343}, {2.09998, 7.25237}, {2.13331,6.95413}, {2.16664, 6.63875}, {2.19998, 6.34642}, {2.23331,6.06828}, {2.26664, 5.77579}, {2.29998, 5.4747}, {2.33331, 5.15976}, {2.36664, 4.84916}, {2.39998, 4.5256}, {2.43331,4.22336}, {2.46664, 3.9177}, {2.49998, 3.58284}, {2.53331,3.2908}, {2.56664, 2.97411}, {2.59998, 2.6861}, {2.63331, 2.4965}, {2.66664, 2.73492}, {2.69998, 2.99366}, {2.73331, 3.29602}, {2.76663, 3.58096}, {2.79997, 3.83507}, {2.8333,4.1179}, {2.86663, 4.39381}, {2.89997, 4.66047}, {2.9333, 4.95059}, {2.96663, 5.23038}, {2.99997, 5.48554}, {3.0333, 5.77507}, {3.06663, 6.03556}, {3.09997, 6.30288}, {3.1333,6.56806}, {3.16663, 6.82612}, {3.19997, 7.11681}, {3.2333,7.37396}, {3.26663, 7.63213}, {3.29997, 7.89755}, {3.3333, 8.15167}, {3.36663, 8.4428}, {3.39997, 8.6969}, {3.4333,8.95516}, {3.46663, 9.22325}, {3.49997, 9.47407}, {3.5333,9.73972}, {3.56663, 9.98549}, {3.59997, 10.2457}, {3.6333,10.4917}, {3.66663, 10.7494}, {3.69997, 10.9985}, {3.7333,11.2493}, {3.76663, 11.5069}, {3.79997, 11.7599}, {3.8333,12.0148}, {3.86663, 12.2645}, {3.89996, 12.5198}, {3.93329,12.7714}, {3.96662, 13.0222}, {3.99996, 13.2753}, {4.03329,13.4973}, {4.06662, 13.7457}, {4.09996, 13.9856}, {4.13329,14.2364}, {4.16662, 14.4828}, {4.19996, 14.7348}, {4.23329,14.9753}, {4.26662, 15.211}, {4.29996, 15.4466}, {4.33329,15.6922}, {4.36662, 15.9198}, {4.39996, 16.1627}, {4.43329,16.4001}, {4.46662, 16.6353}, {4.49996, 16.8629}, {4.53329,17.1011}, {4.56662, 17.3418}, {4.59996, 17.5674}, {4.63329,17.81}, {4.66662, 18.0313}, {4.69996, 18.2533}, {4.73329,18.4823}, {4.76662, 18.7227}, {4.79996, 18.9488}, {4.83329,19.1835}, {4.86662, 19.4019}, {4.89996, 19.6282}, {4.93329,19.86}, {4.96662, 20.084}, {4.99994, 20.3083}, {5.03328,20.5353}, {5.06661, 20.7602}, {5.09994, 20.9745}, {5.13328, 21.1844}, {5.16661, 21.4296}, {5.19994, 21.6461}, {5.23328,21.8579}, {5.26661, 22.0885}, {5.29994, 22.3081}, {5.33328,22.5211}}];

The measurement shows, where/when the collision takes place

{tc, xc} = MinimalBy[tx, Last][[1]];
(*{2.63331, 0.024965}*)

The collision (which isn't measured!) is described by the restitution coefficient x'[SuperPlus[tc]]==-e x'[ SuperMinus[tc]]

Modified system (only describes the state before/after the collision) x''[t] == -F - km x[t] - cm*x'[t] can be solved piecewise

(*before collision*)
X0 = ParametricNDSolveValue[{ x''[t] == -F - km x[t]   - cm*x'[t] , 
x'[tc] == v0 , x[tc] == xc}, x, {t, tx[[1, 1]], tc}, { v0, F, km, cm , e }]

(*after collision*)
X1 = ParametricNDSolveValue[{ x''[t] == -F - km x[t]   - cm*x'[t] , 
x'[tc] == -v0 e, x[tc] == xc}, x, {t, tc, tx[[-1, 1]]}, { v0, F, km, cm, e  }]

system identification

mod=NonlinearModelFit[tx, {Which[t <= tc, X0[v0, F, km, cm , e ][t],t > tc, X1[v0, F, km, cm , e ][t]], 0 < e < 1, F > 0, km > 0,cm > 0}, 
{v0, F, km, cm , e}, t, Method -> "NMinimize"]

shows

Show[{ListPlot[tx, PlotStyle -> Red],Plot[mod[t], {t, 0, tx[[-1, 1]]}]}]

enter image description here

very good agreement with the measurement and justfies the use of a different modell.

$\endgroup$
2
  • $\begingroup$ (+1) Yeah, pretty straightforward, after seeing this. :) $\endgroup$ – Anton Antonov Dec 11 '20 at 15:23
  • $\begingroup$ @AntonAntonov Thanks! But I'm unsure wether it makes sense to additional modell the collision in more detail. $\endgroup$ – Ulrich Neumann Dec 11 '20 at 15:33
2
$\begingroup$
  • This answer does not take into account all details about units and modeled process given by OP.

    • Hence it should be seen as "in principle" answer.
  • It seems that:

    • Further descriptions of the process and model are needed

    • Multiple modifications of the model and its coding have to be made

  • Please see the comments to the question and this answer.


Here is the measured data:

lsData = {{0., 23.6724}, {0.0333333, 23.4316}, {0.0666667, 23.2125}, {0.1, 22.9737}, {0.133333, 22.7191}, {0.166667, 22.4796}, {0.2, 22.2635}, {0.233333, 22.0175}, {0.266667, 21.7774}, {0.3, 21.5224}, {0.333333, 21.3139}, {0.366667, 21.064}, {0.4, 20.8183}, {0.433333, 20.5699}, {0.466667, 20.3129}, {0.5, 20.0644}, {0.533333, 19.8333}, {0.566656, 19.5862}, {0.599989, 19.3391}, {0.633322, 19.094}, {0.666656, 18.8495}, {0.699989, 18.5973}, {0.733322, 18.3451}, {0.766656, 18.09}, {0.799989, 17.8299}, {0.833322, 17.581}, {0.866656, 17.3204}, {0.899989, 17.0659}, {0.933322, 16.817}, {0.966656, 16.5627}, {0.999989, 16.3046}, {1.03332, 16.0535}, {1.06666, 15.7956}, {1.09999, 15.5383}, {1.13332, 15.2806}, {1.16666, 15.0236}, {1.19999, 14.7635}, {1.23332, 14.5015}, {1.26666, 14.2514}, {1.29999, 13.9673}, {1.33332, 13.6998}, {1.36666, 13.4402}, {1.39999, 13.1574}, {1.43332, 12.8848}, {1.46666, 12.6188}, {1.49999, 12.3376}, {1.53332, 12.0596}, {1.56666, 11.7867}, {1.59999, 11.5302}, {1.63332, 11.2418}, {1.66664, 10.9721}, {1.69998, 10.7005}, {1.73331, 10.399}, {1.76664, 10.1111}, {1.79998, 9.83385}, {1.83331, 9.56173}, {1.86664, 9.25114}, {1.89998, 8.98928}, {1.93331, 8.70041}, {1.96664, 8.41822}, {1.99998, 8.13319}, {2.03331, 7.84509}, {2.06664, 7.53343}, {2.09998, 7.25237}, {2.13331, 6.95413}, {2.16664, 6.63875}, {2.19998, 6.34642}, {2.23331, 6.06828}, {2.26664, 5.77579}, {2.29998, 5.4747}, {2.33331, 5.15976}, {2.36664, 4.84916}, {2.39998, 4.5256}, {2.43331, 4.22336}, {2.46664, 3.9177}, {2.49998, 3.58284}, {2.53331, 3.2908}, {2.56664, 2.97411}, {2.59998, 2.6861}, {2.63331, 2.4965}, {2.66664, 2.73492}, {2.69998, 2.99366}, {2.73331, 3.29602}, {2.76663, 3.58096}, {2.79997, 3.83507}, {2.8333, 4.1179}, {2.86663, 4.39381}, {2.89997, 4.66047}, {2.9333, 4.95059}, {2.96663, 5.23038}, {2.99997, 5.48554}, {3.0333, 5.77507}, {3.06663, 6.03556}, {3.09997, 6.30288}, {3.1333, 6.56806}, {3.16663, 6.82612}, {3.19997, 7.11681}, {3.2333, 7.37396}, {3.26663, 7.63213}, {3.29997, 7.89755}, {3.3333, 8.15167}, {3.36663, 8.4428}, {3.39997, 8.6969}, {3.4333, 8.95516}, {3.46663, 9.22325}, {3.49997, 9.47407}, {3.5333, 9.73972}, {3.56663, 9.98549}, {3.59997, 10.2457}, {3.6333, 10.4917}, {3.66663, 10.7494}, {3.69997, 10.9985}, {3.7333, 11.2493}, {3.76663, 11.5069}, {3.79997, 11.7599}, {3.8333, 12.0148}, {3.86663, 12.2645}, {3.89996, 12.5198}, {3.93329, 12.7714}, {3.96662, 13.0222}, {3.99996, 13.2753}, {4.03329, 13.4973}, {4.06662, 13.7457}, {4.09996, 13.9856}, {4.13329, 14.2364}, {4.16662, 14.4828}, {4.19996, 14.7348}, {4.23329, 14.9753}, {4.26662, 15.211}, {4.29996, 15.4466}, {4.33329, 15.6922}, {4.36662, 15.9198}, {4.39996, 16.1627}, {4.43329, 16.4001}, {4.46662, 16.6353}, {4.49996, 16.8629}, {4.53329, 17.1011}, {4.56662, 17.3418}, {4.59996, 17.5674}, {4.63329, 17.81}, {4.66662, 18.0313}, {4.69996, 18.2533}, {4.73329, 18.4823}, {4.76662, 18.7227}, {4.79996, 18.9488}, {4.83329, 19.1835}, {4.86662, 19.4019}, {4.89996, 19.6282}, {4.93329, 19.86}, {4.96662, 20.084}, {4.99994, 20.3083}, {5.03328, 20.5353}, {5.06661, 20.7602}, {5.09994, 20.9745}, {5.13328, 21.1844}, {5.16661, 21.4296}, {5.19994, 21.6461}, {5.23328, 21.8579}, {5.26661, 22.0885}, {5.29994, 22.3081}, {5.33328, 22.5211}};

Below the ODE model programming is changed in several ways:

  • Using RealAbs for x[t]

  • Adding WhenEvent for dealing with the bouncing

  • Using the first x-value of the measurements data to make an initial condition

  • Using parametric formulation for the family of solutions parameterized with k and c

ClearAll[g, m, k, c];
m = 0.3;
g = 9.81;
sol = 
  ParametricNDSolve[{
    m*x''[t] == -k*RealAbs[x[t]]^(3/2) - c*x'[t] - g*m, 
    WhenEvent[x[t] == 0, x'[t] -> -2/3 x'[t]], 
    x'[0] == -3, 
    x[0] == lsData[[1, 2]] 
   }, x, {t, Min[lsData[[All, 1]]], Max[lsData[[All, 1]]]}, {k, c}]

enter image description here

Remark:

  • [...] all we know is that x'[0]==-3, where -3 is the impact velocity before the collision, and x'[T]==2 where 2 is the rebounding velocity after the collision and T is the time of contact, [...]

  • WhenEvent[x[t] == 0, x'[t] -> -2/3 x'[t]] says that when the object touches the ground it bounces (with opposite in sign) velocity that is $2/3$-rds of the velocity just before impact. (The $2/3$ coefficient comes from the velocities described in the question.)


Here we define a function ParDist that measures the deviation of the fit (that takes as arguments parametric function, parameters list, measured data):

Clear[ParDist]
ParDist[x_ParametricFunction, {k_?NumberQ, c_?NumberQ}, tsPath : {{_?NumberQ, _?NumberQ} ..}] := 
   Block[{points, tMin, tMax}, 
    points = Map[{#, x[k, c][#]} &, tsPath[[All, 1]]]; 
    Norm[(tsPath[[All, 2]] - Re[points[[All, 2]]])/tsPath[[All, 2]]] 
   ];

Minimize the measure function ParDist over an appropriate domain for the parameters:

AbsoluteTiming[
  nsol = NMinimize[{ParDist[x /. sol, {k, c}, lsData], -1 <= k <= 0, -2 <= c <= 0}, {k, c}, Method -> "NelderMead", PrecisionGoal -> 3, AccuracyGoal -> 3, MaxIterations -> 100] 
 ]

(* Messages... *)

(*{0.319493, {2.57776, {k -> -0.0223514, c -> -0.0730673}}}*)

(Several experiments can/should be done with different parameter ranges.)


Evaluate the parametric function with the found parameters over the domain of the measured data and plot:

Block[{k, c}, 
   {k, c} = {k, c} /. nsol[[2]]; 
   fitData = Table[{t, Re[x[k, c][t] /. sol]}, {t, lsData[[All, 1]]}] 
  ];
ListPlot[{lsData, fitData}, PlotRange -> All, PlotTheme -> "Detailed",PlotLegends -> {"Measured", "Fitted"}]

enter image description here


Similar, but more complicated procedure is described in this answer of "Model calibration with phase space data".

$\endgroup$
10
  • $\begingroup$ Sorry I think you might have misunderstood my question. The data that I provided was not very useful as we would want to model the instant of collision instead of the entire bouncing motion of the object. $\endgroup$ – bob the legend Dec 10 '20 at 12:08
  • $\begingroup$ My ODE is basically saying that when the object barely touches the ground with some initial velocity v, x=0. The object will then deform in shape, thus "sinking" into the ground where x<0. We are trying to model this deformation depth using the equation. All we know is the velocity at the start, and the velocity at the end. We know roughly what the time of contact is but still unsure. $\endgroup$ – bob the legend Dec 10 '20 at 12:11
  • $\begingroup$ But thank you anyways because the answer is still very useful as you have taught me the basic method for fitting the data :) $\endgroup$ – bob the legend Dec 10 '20 at 12:12
  • 1
    $\begingroup$ @AntonAntonov It is not clear, how you correlated data for x given in cm with equation where free fall acceleration is given in SI unit $g=9.81 m/s^2$? :) $\endgroup$ – Alex Trounev Dec 10 '20 at 12:42
  • 1
    $\begingroup$ @AntonAntonov It's ok and (+1) for optimization. Actually we can neglect by g in this problem or use it as parameter to be optimized too. $\endgroup$ – Alex Trounev Dec 10 '20 at 13:59
2
$\begingroup$

This is an extension for the excellent answer from @Ulrich Neumann considering

$$m\ddot x=-kx^{\alpha}-c\dot x-mg$$ instead of

$$m\ddot x=-kx-c\dot x-mg$$

tx = Map[{#[[1]], #[[2]]/100} &, data]
{tc, xc} = MinimalBy[tx, Last][[1]];

X0 = ParametricNDSolveValue[{x''[t] == -F - km Sign[x[t]] Abs[x[t]]^alpha - cm*x'[t], x'[tc] == v0, x[tc] == xc}, x, {t, tx[[1, 1]], tc}, {v0, F, km, cm, alpha, e}]
X1 = ParametricNDSolveValue[{x''[t] == -F - km Sign[x[t]] Abs[x[t]]^alpha - cm*x'[t], x'[tc] == -v0 e, x[tc] == xc}, x, {t, tc, tx[[-1, 1]]}, {v0, F, km, cm, alpha, e}]

mod = NonlinearModelFit[tx, {Which[t <= tc, X0[v0, F, km, cm, alpha, e][t], t > tc, X1[v0, F, km, cm, alpha, e][t]], 0 < e < 1, F > 0, km > 0, cm > 0, 0.5 < alpha < 3}, {v0, F, km, cm, alpha, e}, t, Method -> "NMinimize"]

Show[{ListPlot[tx, PlotStyle -> Red], Plot[mod[t], {t, 0, tx[[-1, 1]]}]}]

Normal[mod]

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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