9
$\begingroup$

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}]
$\endgroup$
1

1 Answer 1

12
$\begingroup$

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

$\endgroup$
8
  • $\begingroup$ 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]]. $\endgroup$ Commented Feb 14, 2018 at 10:01
  • $\begingroup$ 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. $\endgroup$
    – joojaa
    Commented Feb 14, 2018 at 10:14
  • $\begingroup$ @joojaa, very good point., thank you. Just took the setup in the docs without thinking about it. $\endgroup$
    – kglr
    Commented Feb 14, 2018 at 10:16
  • 1
    $\begingroup$ @SjoerdSmit, you are right. ParametricNDSolve my first try; but somehow i couldn't get it right. I went with the example in the docs. $\endgroup$
    – kglr
    Commented Feb 14, 2018 at 10:22
  • 1
    $\begingroup$ 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 $\endgroup$
    – joojaa
    Commented Feb 14, 2018 at 10:45

Your Answer

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

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