1
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ You could start by trying to fit a simpler equation. There are a lot of errors there $\endgroup$ – Dr. belisarius Dec 12 '14 at 18:38
  • $\begingroup$ @belisarius Do you mean syntax errors? $\endgroup$ – mattiav27 Dec 12 '14 at 19:22
  • $\begingroup$ No. I mean conceptual errors. Like your orbita[] definition and usage. You need more training on easier things $\endgroup$ – Dr. belisarius Dec 12 '14 at 19:24
  • $\begingroup$ 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. $\endgroup$ – Michael E2 Dec 12 '14 at 19:25
  • $\begingroup$ @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 $\endgroup$ – Dr. belisarius Dec 12 '14 at 19:32
3
$\begingroup$

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]]

Mathematica graphics

It does not look too bad.

$\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.