4
$\begingroup$

I have written the following code to fit the function g2 to data (a,b,c,d) but I don't realize where I have made a mistake,without d is working well but when I put parameter d in exponent it does not work.

data = {{6.47, 3.65}, {7.43, 
    3.45}, {3.9, -2.94}, {4.8, -1.29}, {2.48, -0.35}, {6.32, 
    3.16}, {2.59, -1.19}, {9.13, -2.}, {3.81, -3.04}, {3.33, -2.68}};  

g2[a_?NumberQ, b_?NumberQ, c_?NumberQ] :=    
 First[y /. NDSolve[{y''[x] + a y[x] == 0, y[0] == b, y'[0] == c}, y, {x, 0, 10}]] 

nlm2 = NonlinearModelFit[data, g2[a, b, c][x]^d, {a, b, c, d}, x]

Could you please advise me why this code does not work? Thank you so much

|improve this question
$\endgroup$
  • 4
    $\begingroup$ What is g2^r[x] supposed to be? Maybe you wanted g2[a, b, c][x]^d? But, consider also the use of ParametricNDSolve[]. $\endgroup$ – J. M. will be back soon♦ Dec 29 '16 at 22:28
  • $\begingroup$ g2 has defined above, r is constant which is suppose to obtain from fitting $\endgroup$ – amin bk Dec 29 '16 at 22:34
  • 3
    $\begingroup$ This is not how it works in Mathematica. If g[x] is a function, then the syntax to square it is g[x]^2 and not g^2[x]. This is what JM was referring to. $\endgroup$ – Nasser Dec 29 '16 at 22:41
  • $\begingroup$ You need to define how d enters the equation. It's currently in your parameter list but nowhere else. If r is another constant to be estimated, it also needs to be in the parameter list and entered into the function correctly (as was mentioned by @J.M.). $\endgroup$ – JimB Dec 29 '16 at 23:45
  • 1
    $\begingroup$ g2 has an explicit result using DSolve (b Cos[Sqrt[a] x] + (c Sin[Sqrt[a] x])/Sqrt[a]) so that should simplify things. However, the recent edit by @aminbk raises that result to a potentially non-integer power (d). If that is correct, doesn't d need to be restricted to integer powers to get a real number? $\endgroup$ – JimB Dec 30 '16 at 0:42
6
$\begingroup$

Maybe so:

data = {{2.48, -0.35}, {2.59, -1.19}, {3.33, -2.68}, {3.81, -3.04},
{3.9, -2.94}, {4.8, -1.29}, {6.32, 3.16}, {6.47, 3.65}, {7.43, 3.45}, {9.13, -2.}};

pfun = ParametricNDSolveValue[{y''[x] + a*y[x] == 0, y[0] == b, y'[0] == c}, y, {x, 0, 10}, {a, b, c}];

fit = FindFit[data, pfun[a, b, c][x], {{a, 0.1}, {b, 0.2}, {c, 0.5}},x]

(*{a -> 1.04037, b -> 2.47101, c -> 2.29589}*)

Plot[pfun[a, b, c][x] /. fit, {x, 0, 10}, Epilog -> {PointSize[Medium], Point@data}]

enter image description here

DSolveValue[{y''[x] + a*y[x] == 0, y[0] == b, y'[0] == c}, y[x], x] /. fit

$0.980404 (2.29589 \sin (1.01999 x)+2.5204 \cos (1.01999 x))$

|improve this answer
$\endgroup$
  • $\begingroup$ How about the parameter "d"? $\endgroup$ – amin bk Dec 29 '16 at 23:13
  • 4
    $\begingroup$ @amin, it's not Mariusz's fault that there was no mention of how d enters the function to be fit against in the question. $\endgroup$ – J. M. will be back soon♦ Dec 29 '16 at 23:52
3
$\begingroup$

There is an explicit solution for the differential equation which can simplify things:

DSolve[{y''[x] + a y[x] == 0, y[0] == b, y'[0] == c}, y, {x, 0, 10}][[1, 1, 2, 2]]

b Cos[Sqrt[a] x] + (c Sin[Sqrt[a] x])/Sqrt[a]

If the recent edit in the original question really requires that function raised to a power that must be estimated with the data, then that value can only take on odd integer values (so that we can obtain real and sometimes negative values to match the data).

NonlinearModelFit doesn't allow such restrictions of parameters (as far as I can tell) but we can set various values of d and choose the value of d that minimizes either the $AIC$ or $AIC_c$ value.

d = {-9, -7, -5, -3, -1, 1, 3, 5, 7, 9}
aicc = Table["aicc", {i, Length[d]}];
Do[aicc[[i]] = 
   NonlinearModelFit[data, {(b Cos[Sqrt[a] x] + (c Sin[Sqrt[a] x])/Sqrt[a])^d[[i]]},
     {a, b, c}, x]["AICc"],
      {i, Length[d]}];

ListPlot[{{{1, aicc[[n + 1]]}}, Transpose[{d, aicc}]}, 
 PlotStyle -> {{Green, PointSize[0.03]}, Red}]

AICc values for d

We see that d = 1 results in the smallest $AIC_c$ value.

So to estimate the remaining parameters we set d = 1 and find the 95% confidence bands for the fit.

nlm = NonlinearModelFit[data, 
   b Cos[Sqrt[a] x] + (c Sin[Sqrt[a] x])/Sqrt[a], {a, b, c}, x];

nlm["BestFitParameters"]
(* {a -> 1.0403738146939832, b -> 2.4710141618669903, c -> 2.295887035294033} *)

mpb = nlm["MeanPredictionBands"];

Show[ListPlot[data],
 Plot[{nlm[x], mpb}, {x, Min[data[[All, 1]]], Max[data[[All, 1]]]}]]

Plot of data, fit, and 95% confidence bands

The residuals aren't great as the points near the peak and the trough are all underestimated. Another model might be considered or an examination into why the data isn't a great fit to the theoretical model. (Note that the confidence bands don't account for the fact that we needed to estimate d.)

|improve this answer
$\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.