How can I fit a differential equation to experimental data?

Question

I would like to add a fit to the following data. Here is the differential equation I am using. g=9.81 l=0.2 m=0.09

{((-g*Sin[y[x]])/(l)) - ((k)(y'[x]))/(m) == (y''[x]), y[0] == 1.22173, y'[0] == 0}, y, {x, 0, 3*Pi}

I need to find values for k.

Here is the data

I've also added the data over here.

You can copy-paste this data into mathematica, and it will work.

@anderstood I believe you could just copy-paste it into mathematica? Does that work? I tried pasting the data here, but slashes appear at the end of each line. — fidafa123, Feb 18 '18 at 13:36
I am trying to find a coefficient for k. g is 9.81, and m is 0.9 — fidafa123, Feb 18 '18 at 13:52
@anderstood Would that equation work for pendulums at large angles? Because the equation I have uses sintheta. — fidafa123, Feb 18 '18 at 13:53
I have added the data, so are you able to copy paste it into mathematica? — fidafa123, Feb 18 '18 at 14:12
Yes, just l is missing now. You should add them in the equation. — anderstood, Feb 18 '18 at 14:13
l=0.2, and g is 9.81. I need mathematica to find a value for k, and graph the fit. — fidafa123, Feb 18 '18 at 14:16
The data can be retrieved with Most@ToExpression[ "{" <> Import["http://m.uploadedit.com/bbtc/1518933165520.txt"] <> "}"] -- there are minor syntax errors that are fixed with the braces and Most, at least currently. — Michael E2, Feb 18 '18 at 14:48
Btw k for damping is a terrible notation... ($k$ usually denotes stiffness). — anderstood, Feb 18 '18 at 15:28
Ahh, its supposed to be b actually, in many books, but I usually use k for just about any constant. — fidafa123, Feb 18 '18 at 15:30

OkkesDulgerci · Answer 1 · 2018-02-19T02:21:37.683

2

When we assume g,l and m are known fit is not good.

g = 9.81; l = 0.2 ; m = 0.09;

model[k_?NumberQ] := (model[k] = 
   Module[{y, x}, 
    NDSolveValue[{y''[x] == -g*Sin[y[x]]/l - k y'[x]/m, 
      y[0] == data[[1, 2]], y'[0] == 0}, y, {x, 0, 3*Pi}]])

nlm = NonlinearModelFit[data, model[k][x], k, x, 
   Method -> {NMinimize, Method -> "DifferentialEvolution"}];

Plot[nlm[x], {x, 0, data[[-1, 1]]}, 
 Epilog -> {Red, PointSize[Medium], Point[data]}, PlotRange -> All, 
 Frame -> True]
nlm["BestFitParameters"]

{k -> 0.405883}

We can do better than this when we assume we only know g.

 g = 9.81;

model[k_?NumberQ, l_?NumberQ, 
  m_?NumberQ] := (model[k, l, m] = 
   Module[{y, x}, 
    NDSolveValue[{y''[x] == -g*Sin[y[x]]/l - k y'[x]/m, 
      y[0] == data[[1, 2]], y'[0] == 0}, y, {x, 0, 3*Pi}]])

nlm = NonlinearModelFit[data, model[k, l, m][x], {k, l, m}, x, 
   Method -> {NMinimize, Method -> "DifferentialEvolution"}];

Plot[nlm[x], {x, 0, data[[-1, 1]]}, 
 Epilog -> {Red, PointSize[Medium], Point[data]}, PlotRange -> All, 
 Frame -> True]
nlm["BestFitParameters"]

{k -> 1.25031, l -> 0.298879, m -> 10.5371}

g = 9.81;

model[k_?NumberQ, 
  l_?NumberQ] := (model[k, l] = 
   Module[{y, x}, 
    NDSolveValue[{y''[x] + k y'[x] == -g*Sin[y[x]]/l, 
      y[0] == data[[1, 2]], y'[0] == 0}, y, {x, 0, 3*Pi}]])

nlm = NonlinearModelFit[data, model[k, l][x], {k, l}, x, 
   Method -> {NMinimize, Method -> "DifferentialEvolution"}];

Plot[nlm[x], {x, 0, data[[-1, 1]]}, 
 Epilog -> {Red, PointSize[Medium], Point[data]}, PlotRange -> All, 
 Frame -> True]
nlm["BestFitParameters"]

{k -> 0.119336, l -> 0.29895}

edited Feb 19 '18 at 02:21

answered Feb 18 '18 at 15:11

OkkesDulgerci

10,716
1
19
38

Hi, I entered your code, but I get the following errors: http://prntscr.com/igfays – fidafa123 Feb 18 '18 at 15:16
Also, does it give a value for k? g, l, m are known constants, with g=9.81, l=0.2 and m=0.09 – fidafa123 Feb 18 '18 at 15:18
l and m are wrong values. So I fit them also. These are the values {k -> 1.25031, l -> 0.298879, m -> 10.5371} – OkkesDulgerci Feb 18 '18 at 15:20
Actually m is very off. l is fine.. I don't know how you find m – OkkesDulgerci Feb 18 '18 at 15:22
Remove l and m from your code and use Clear[l,m] – OkkesDulgerci Feb 18 '18 at 15:24
I see, what code did you use to find k? – fidafa123 Feb 18 '18 at 15:24
I think that is due to the fact that my experiment may not have been carried out well. It is to be noted that k/m is supposed to be a very small value, since there is supposed to be little to no friction. – fidafa123 Feb 18 '18 at 15:27
maybe instead, remove the 'm' altogether? Then finally divide the k value by m? – fidafa123 Feb 18 '18 at 15:28
m, g and l are supposed to be constants – fidafa123 Feb 18 '18 at 15:29
I actually did my experiment using a length of 0.2meters, and a mass of 0.09 kilos. Would it be possible to remove the m, and the make the entire coefficient in front of y'[x] equal to k? – fidafa123 Feb 19 '18 at 01:05
I did it manually, and let m=0.3 instead. https://prnt.sc/igmn2j Would it be possible to make mathematica find the best value of k? – fidafa123 Feb 19 '18 at 01:07
You mean l=0.3? – OkkesDulgerci Feb 19 '18 at 01:49
sorry, yes, l=0.3 – fidafa123 Feb 19 '18 at 03:16
I meant as in g and l are known, and the code is supposed to find a k-value. – fidafa123 Feb 19 '18 at 03:18
@OmkarVaidya this is probably coming too late to help you. But just glancing at the results vs your known values you probably have a unit conversion error. I'm assuming m is mass in this case. So if it's off by a factor of 1000, that implies perhaps a gram vs kilogram unit mixup. – OmnipotentEntity Jan 24 '19 at 21:27

score 1 · Answer 2 · answered Feb 18 '18 at 14:44

The following finds the best fit (in some sense defined in the cost function). However based on the data and model you provided, the result is poor --- which is good for you, because having the solution directly from someone else is not very interesting anyway...

l = 0.2; g = 9.81; m = 0.9;
eq[k_] := -g/l*Sin[y[x]] - k*y'[x]/m == y''[x]; 
inter = Interpolation[data];
ic = {y[0] == inter[0], y'[0] == inter'[0]};
cost[k_?NumericQ] := Block[{sol},
  sol = NDSolveValue[{eq[k]}~Join~ic, y, {x, 0, 3}];
  NIntegrate[(sol[x] - inter[x])^2, {x, 0., 3}]]
FindMinimum[cost[k], {k, 4}]

is it possible to graph this over the data? – fidafa123 Feb 18 '18 at 14:48 — fidafa123, Feb 18 '18 at 14:48

How can I fit a differential equation to experimental data?

2 Answers2

Linked