9

I've done an experiment where I swung a pendulum under air resistance. Is it possible to model the data using the following differential equation and find a b-value?

(y''[x])+ Sin[y[x]] + b(y'[x]) == 0, y[0] == 1.5, y'[0] == 0},  y, {x, 0, 3*Pi}]
m_goldberg
  • 107,779
  • 16
  • 103
  • 257
fidafa123
  • 145
  • 1
  • 8

1 Answers1

12

Mimicking the examples in 

ClearAll[x, y, b, β, model]
b0 = .7;
sol = First[y /. NDSolve[{y''[x] + Sin[y[x]] + b0  y'[x] == 0, y[0] == 1.5, y'[0] == 0}, 
  y, {x, 0, 3 Pi}]];
xvals = N[Range[0, 3 Pi, 3 Pi/100]];
data = Transpose[{xvals, sol[xvals] + RandomReal[{-.1, .1}, 101]}];

enter image description here

FindFit

model[b_?NumberQ] := (model[b] = First[y /. 
 NDSolve[{y''[x] + Sin[y[x]] + b (y'[x]) == 0, y[0] == 1.5, y'[0]==0}, y, {x, 0, 3 Pi}]])

fit = FindFit[data, model[β][x], {{β, .1}}, x, PrecisionGoal -> 4, AccuracyGoal -> 4]

{β -> 0.695487}

Show[ListPlot[data], Plot[model[β][x] /. fit, {x, 0, 3 Pi}, PlotStyle -> Green]]

enter image description here

NonlinearModelFit

nlm = NonlinearModelFit[data, model[β][x], {{β, .1}}, x, 
  PrecisionGoal -> 4, AccuracyGoal -> 4];
Show[ListPlot[data], Plot[nlm[x], {x, 0, 3 Pi}, PlotStyle -> Red]]

enter image description here

An alternative (4-parameter) model:

ClearAll[model]
model[a_?NumberQ, b_?NumberQ, c_?NumberQ, d_?NumberQ] := (model[a, b, c, d] = 
  First[y /. NDSolve[{y''[x] + a Sin[b  y[x]] + c (y'[x]) == 0, y[0] == d, y'[0] == 0}, 
   y, {x, 0, 3 Pi}]])
nlm = NonlinearModelFit[data, model[α, β, γ, δ ][x], 
 {{α, .1}, {β, .1}, {γ, .1}, {δ, .1}}, x, PrecisionGoal -> 4, AccuracyGoal -> 4];
nlm["ParameterTable"]

enter image description here

Show[ListPlot[data], Plot[nlm[x], {x, 0, 3 Pi}, PlotStyle -> Purple]]

enter image description here

kglr
  • 394,356
  • 18
  • 477
  • 896
  • ParametricNDSolve can come in handy too. For example, the first fit will also work with model = ParametricNDSolve[{y''[x] + Sin[y[x]] + b (y'[x]) == 0, y[0] == 1.5, y'[0] == 0}, y, {x, 0, 3 Pi}, b][[1, 2]]. – Sjoerd Smit Feb 14 '18 at 10:01
  • The plots would probably be nicer if your error would be symmetric about the data (-0.05 to 0.05 instead of 0 to 0.1) and not just on the positive side of the data. – joojaa Feb 14 '18 at 10:14
  • @joojaa, very good point., thank you. Just took the setup in the docs without thinking about it. – kglr Feb 14 '18 at 10:16
  • 1
    @SjoerdSmit, you are right. ParametricNDSolve my first try; but somehow i couldn't get it right. I went with the example in the docs. – kglr Feb 14 '18 at 10:22
  • 1
    And now the fits are near perfect too. Lesson to be learned when fitting real data allways put a + offset into the fit to find any possible systematic error – joojaa Feb 14 '18 at 10:45
  • @kglr I pasted the code into mathematica. However, I get the following graph. I also changed the initial condition of the function to be 1.22173. https://prnt.sc/igaudf – fidafa123 Feb 18 '18 at 05:47
  • @kglr Here is my data: http://m.uploadedit.com/bbtc/1518933165520.txt – fidafa123 Feb 18 '18 at 05:54
  • @kglr its this link actually: http://prntscr.com/igaxb8 – fidafa123 Feb 18 '18 at 06:02