# FindFit::nrlnum: error with FindFit

I am trying to write a code to fit actual data of the orbit of a star with Mathematica 8.

First I tried to write a code to fit a mock orbit, using the Application with differential equations (the first one) as a guide (see here): in that case there is only 1 function I need two (x,y positions) for my orbit.

Here is my attempt up to now:

Msun = 2*10^30;
M = 4.3*10^6;
G = 6.67*10^-11;
UA = 150*10^6;
(*mock orbit*)
s = NDSolve[{Vx'[t] == (-M*G*X[t]*Msun)/(X[t]^2 + Y[t]^2)^(3/2),
Vy'[t] == (-M*G*Y[t]*Msun)/(X[t]^2 + Y[t]^2)^(3/2),
X'[t] == Vx[t], Y'[t] == Vy[t], Vx == #1, Vy == #2,
X == #3, Y == #4}, {Vx, Vy, X, Y}, {t, 0, 10^4}] &;

tempi = Range[0, 10^4, 10^2];
Xdat = Flatten[
Evaluate[X[#] /. s[0, 4*10^7, 1000*UA, 500*UA]] & /@ tempi, 1];
Ydat = Flatten[
Evaluate[Y[#] /. s[0, 4*10^7, 1000*UA, 500*UA]] & /@ tempi, 1];
dati = Transpose[{Xdat*, Ydat}];
ListPlot[dati]

(*this is the model*)
orbita[a_?NumberQ, b_?NumberQ, c_?NumberQ, d_?NumberQ,
e_?NumberQ] := (orbita[a, b, c, d, e] := {X, Y} /.
NDSolve[{Vx'[t] == (-e*G*X[t]*Msun)/(X[t]^2 + Y[t]^2)^(3/2),
Vy'[t] == (-e*G*Y[t]*Msun)/(X[t]^2 + Y[t]^2)^(3/2),
X'[t] == Vx[t], Y'[t] == Vy[t], Vx == a, Vy == b,
X == c, Y == d}, {Vx, Vy, X, Y}, {t, 0, 10^4}]);

FindFit[dati,
orbita[0, 4*10^7, 1000*UA, 500*UA, Mc][X, Y], {{Mc, 10^6}}, t]


I get the following error:

FindFit::nrlnum:"The function value {-7.97408*10^10+Null[{X,Y}],<<49>>,<<51>>} is not \ a list of real numbers with dimensions {101} at {Mc} = {1.*^6}."

I have also tried with:

FindFit[dati,
orbita[0, 4*10^7, 1000*UA, 500*UA, Mc][X, Y], {{Mc, 10^6}}, {X,Y}]


but I get another error:

FindFit::fitc: Number of coordinates (1) is not equal to the number of variables (2)

Where are my errors?

• You could start by trying to fit a simpler equation. There are a lot of errors there Dec 12, 2014 at 18:38
• @belisarius Do you mean syntax errors? Dec 12, 2014 at 19:22
• No. I mean conceptual errors. Like your orbita[] definition and usage. You need more training on easier things Dec 12, 2014 at 19:24
• You have := instead of = in second eq. in def. of orbita; you should have {X[t], Y[t]} /. NDSolve... and no [X, Y] in FindFit. But the real question is whether FindFit can fit 2D output data as a function of 1D input. It seems to do only multivariate 1D fitting; you want univariate 2D. Dec 12, 2014 at 19:25
• @MichaelE2 But for example Mc, the fitted var isn't used anywhere, and a few more errors too. The real question could be how you should approach a complex problem: one step at a time Dec 12, 2014 at 19:32

First, some slight changes to orbita:

orbita[a_?NumberQ, b_?NumberQ, c_?NumberQ, d_?NumberQ,
e_?NumberQ] := (orbita[a, b, c, d, e] = {X[t], Y[t]} /.
First@NDSolve[{Vx'[t] == (-e*G*X[t]*Msun)/(X[t]^2 + Y[t]^2)^(3/2),
Vy'[t] == (-e*G*Y[t]*Msun)/(X[t]^2 + Y[t]^2)^(3/2),
X'[t] == Vx[t], Y'[t] == Vy[t], Vx == a, Vy == b,
X == c, Y == d}, {Vx, Vy, X, Y}, {t, 0, 10^4}]);


Memoization (or caching; see FindFit documentation: 1 2) is done with =; the := causes the return value to be Null and is the source of your first error message. The other problem is the what to return. I suggest {X[t], Y[t]}, with the variable t in place. I'd also strip an extra set of braces {} with First@NDSolve....

Second, I don't think FindFit will work on 2D univariate data ({x, y} as a function of t). At least I could find no example and a naive toy trial failed. So use FindMinimum to minimize the sum of squares. The objective function is given by

ClearAll[obj];
obj[Mc_?NumericQ] :=
Total[#.# & /@ ((orbita[0, 4*10^7, 1000*UA, 500*UA, Mc] /. t -> # & /@
tempi) - dati)]


Then minimize:

{min, sol} = FindMinimum[obj[Mc], {Mc, 10^6}]


FindMinimum::lstol: The line search decreased the step size to within the tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the function. You may need more than MachinePrecision digits of working precision to meet these tolerances. >>

(*  {1.54869*10^9, {Mc -> 4.3*10^6}}  *)


Inspect the solution to see if the warning is significant.

Show[
ParametricPlot[
orbita[0, 4*10^7, 1000*UA, 500*UA, Mc] /. sol, {t, 0, Max[tempi]},
PlotStyle -> None, Mesh -> {tempi},
MeshStyle -> {PointSize[Large], Red}],
ListPlot[dati]]
` It does not look too bad.