Finding a good surface fitting function
- To: mathgroup at smc.vnet.net
- Subject: [mg116819] Finding a good surface fitting function
- From: Marco Masi <marco.masi at ymail.com>
- Date: Tue, 1 Mar 2011 05:23:12 -0500 (EST)
I'm trying to find a decent fitting function of a z=f(x,y) surface (ideally a polinomial). I tried with different types of linear combinations of fit functions but get always large deviations between the 'real' surface and the fitting surface Mathematica suggests. Please copy and past the appended code below (runs only on Mathematica 8, and since it evaluates double integrals it is a bit long in execution, be patient) and look at the graphical output which shows the differences. The important line is that which uses "Fit" (maybe some other works better?) Can anyone achieve better fits (ideally less than 1% overall deviation)?
Marco.
------------------------------------------------------------
(* Fundamental parameter setting *)
pc = 3600*24*365*299792458*3.2638; M0 = 1.989 10^30; Myr =
3600*24*365*10^6;
\[Alpha] = (M0*Myr^2)/pc^3; G = \[Alpha] * 6.672 10^-11;
\[CapitalSigma] = 492; rd = 3500; (* Central density and scale radius \
of the Freeman disk *)
Subscript[\[Phi], Free][X_, Y_] := - \[Pi] G \[CapitalSigma] Sqrt[
X^2 + Y^2] (BesselI[0, Sqrt[X^2 + Y^2]/(2 rd)] BesselK[1, Sqrt[
X^2 + Y^2]/(2 rd)] -
BesselI[1, Sqrt[X^2 + Y^2]/(2 rd)] BesselK[0, Sqrt[X^2 + Y^2]/(
2 rd)]);
Subscript[\[Phi], Free3D][R_,
Z_] := - 2 \[Pi] G \[CapitalSigma] rd^2 NIntegrate[
Exp[-k Abs[Z]]/(1 + (k rd)^2)^1.5 (BesselJ[0, k R]), {k,
0, \[Infinity]}, Method -> "LevinRule"];
Rmin = -15001; Rmax = 15000; ; Rstep = 1000;
Zmin = -15000; Zmax = 15000; Zstep = 1000;
PotRange = -100000 ;
Print["Building the FreemanTable. Wait..."];
FreemanTable =
Parallelize[
Table[Subscript[\[Phi], Free3D][R, Z], {R, Rmin, Rmax, Rstep}, {Z,
Zmin, Zmax, Zstep}]];
RDim = Part[Dimensions[FreemanTable], 1]; ZDim =
Part[Dimensions[FreemanTable], 2];
GalPotPoints =
Flatten[Table[{Rmin + Rstep*(R - 1), Zmin + Zstep*(Z - 1),
FreemanTable[[R, Z]]}, {R, 1, RDim}, {Z, 1, ZDim}], 1];
MinPot = Min[
Table[GalPotPoints[[i, 3]], {i, 1, Length[GalPotPoints]}]];
ListPointPlot3D[GalPotPoints, PlotRange -> {MinPot, 0},
PlotStyle -> RGBColor[0, 0, 1],
AxesLabel -> {"R", "Z", "\[CapitalPhi]"}, ImageSize -> 800]
Print["The above figure: plot of the Galactic potential points mesh. \
Wait..."];
quad = Fit[
GalPotPoints, {MinPot, X, Y, X*Y, X^2, Y^2, X^3, Y^3, X^2*Y, X*Y^2,
X^4, Y^4, X^3*Y, X^2*Y^2, X*Y^3, Log[Sqrt[X^2 + Y^2]]}, {X, Y}];
Show[Plot3D[quad, {X, Rmin, Rmax}, {Y, Zmin, Zmax},
PlotStyle -> Directive[PointSize[Medium], Blue],
PlotRange -> {MinPot, 0}, AxesLabel -> {"R", "Z", "\[CapitalPhi]"},
ImageSize -> 800],
Graphics3D[{Red, PointSize[0.005], Map[Point, GalPotPoints]}]]
Print["The above figure: plot of the Galactic potential fit \
function", quad, ". Wait..."];
DiffFitGalPotPoints =
Table[100*(1 - (quad /. X -> GalPotPoints[[nr, 1]] /.
Y -> GalPotPoints[[nr, 2]])/GalPotPoints[[nr, 3]]), {nr, 1,
Length[GalPotPoints]}];
DiffFitGalPotPointsPlot =
Table[{GalPotPoints[[nr, 1]], GalPotPoints[[nr, 2]],
DiffFitGalPotPoints[[nr]]}, {nr, 1, Length[GalPotPoints]}];
ListPointPlot3D[DiffFitGalPotPointsPlot,
PlotRange -> {MinDiffPot, MaxDiffPot},
PlotStyle -> RGBColor[1, 0, 1], AxesLabel -> {"R", "Z", "Error %"},
ImageSize -> 800]
Print["The above figure: plot of the difference between the Galactic \
potential and the Galactic potential fit function."];