FindFit::nrlnum: error with FindFit

Question

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[0] == #1, Vy[0] == #2, 
     X[0] == #3, Y[0] == #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[0] == a, Vy[0] == b, 
       X[0] == c, Y[0] == 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?

Answer 1

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[0] == a, Vy[0] == b, 
        X[0] == c, Y[0] == 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.