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}]
$\begingroup$
$\endgroup$
1
-
$\begingroup$ duplicate?: How to fit 3 data sets to a model of 4 differential equations? $\endgroup$– Kuba ♦Feb 14, 2018 at 7:04
Add a comment
|
1 Answer
$\begingroup$
$\endgroup$
8
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]}];
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]]
NonlinearModelFit
nlm = NonlinearModelFit[data, model[β][x], {{β, .1}}, x,
PrecisionGoal -> 4, AccuracyGoal -> 4];
Show[ListPlot[data], Plot[nlm[x], {x, 0, 3 Pi}, PlotStyle -> Red]]
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"]
Show[ListPlot[data], Plot[nlm[x], {x, 0, 3 Pi}, PlotStyle -> Purple]]
-
$\begingroup$
ParametricNDSolvecan come in handy too. For example, the first fit will also work withmodel = 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$ 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$– joojaaFeb 14, 2018 at 10:14
-
$\begingroup$ @joojaa, very good point., thank you. Just took the setup in the docs without thinking about it. $\endgroup$– kglrFeb 14, 2018 at 10:16
-
1$\begingroup$ @SjoerdSmit, you are right.
ParametricNDSolvemy first try; but somehow i couldn't get it right. I went with the example in the docs. $\endgroup$– kglrFeb 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$– joojaaFeb 14, 2018 at 10:45