Fitting complicate functions
- To: mathgroup at smc.vnet.net
- Subject: [mg111049] Fitting complicate functions
- From: Marco Masi <marco.masi at ymail.com>
- Date: Sun, 18 Jul 2010 01:04:59 -0400 (EDT)
I have a very complicate function which has to be inserted into a loop that calculates it zillions of times. This slows the simulation speed and I'm wondering if it is possible to fit the function with a simpler one by adjusting its fitting parameters. I used "Fit" of Mathematica and was partially succesful: however, I tried several functions (polinomial, logarithimic, etc.) and parameters but the difference between the original and the fitting function always tends to diverge in some regions. I'm wondering if this can be done better? The program is appended, let me know if you have a better idea and/or paramterizing fitting function or method.
Regards, Mark.
----------------------------------------------------------------------
(*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; Rmax = 100000;
Subscript[M, BH] = 10^6;
Subscript[M, BG] = 1.6*10^10; rc = 420;
\[CapitalSigma] = 492; rd = 3500;(*Central density and scale radius \
of the Freeman disk*)\[Rho] = 0.05; rh = 4600;(*Central density and \
scale radius of the Psiso halo*)\[Rho]mpsiso = 0.014; rmpsiso = \
13000;(*Central density and scale radius of MPSISO*)\[Rho]nfw = \
0.012; rnfw = 15000;(*Central density and scale radius of NFW*)(*The \
Miyamoto-Nagai potential parameters.*)M1 = 2.4 10^9; M2 =
18.7 10^9; M3 = 40.2 10^9; M4 = 3.4 10^9;
a1 = 80; a2 = 500; a3 = 4600; a4 = 9000; b1 = 95; b2 = 200; b3 =
= 400; \
b4 = 200;
Subscript[\[Phi], BH][X_,
Y_] := -((G Subscript[M, BH])/
Sqrt[X^2 + Y^2]);(*Potential of the Black Hole*)
Subscript[\[Phi], BG][X_,
Y_] := -((G Subscript[M, BG])/Sqrt[X^2 + Y^2 + rc^2]);
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], PSISOMOD][X_,
Y_] := -4 \[Pi] G \[Rho]mpsiso rmpsiso^2 (rmpsiso/
Sqrt[X^2 + Y^2] ArcSinh[Sqrt[X^2 + Y^2]/rmpsiso] -
1/Sqrt[1 + Rmax^2/rmpsiso^2]);
(*Subscript[\[Phi],tot][X_,Y_]=Subscript[\[Phi],BH][X,Y]+Subscript[\
\[Phi],BG][X,Y]+Subscript[\[Phi],Free][X,Y]+Subscript[\[Phi],PSISOMOD]\
[X,Y];*)
Subscript[\[Phi], tot][X_, Y_] = Subscript[\[Phi], Free][X, Y];
Rrange = 15000; PotRange = -50000; Rstep = 301;
Plot3D[Subscript[\[Phi], tot][X, Y], {X, -Rrange,
Rrange}, {Y, -Rrange, Rrange}, PlotRange -> {PotRange, 0},
PlotStyle -> RGBColor[1, 0, 0]]
Print["The above figure: plot of the Galactic potential. Wait..."];
GalPotPoints =
Flatten[Table[{X, Y, Subscript[\[Phi], tot][X, Y]}, {X, -Rrange,
Rrange, Rstep}, {Y, -Rrange, Rrange, Rstep}], 1];
MinPot = Min[
Table[GalPotPoints[[i, 3]], {i, 1, Length[GalPotPoints]}]];
ListPointPlot3D[GalPotPoints, PlotRange -> {MinPot, 0},
PlotStyle -> RGBColor[0, 0, 1]]
Print["The above figure: plot of the Galactic potential points mesh. \
Wait..."];
(* Here "quad" is the fitting function *)
quad = Fit[GalPotPoints, {1, Log[Sqrt[X^2 + Y^2]]}, {X, Y}];
Show[Plot3D[quad, {X, -Rrange, Rrange}, {Y, -Rrange, Rrange},
PlotStyle -> Directive[PointSize[Medium], Blue],
PlotRange -> {MinPot, 0}],
Graphics3D[{Red, PointSize[0.005], Map[Point, GalPotPoints]}]]
Print["The above figure: plot of the Galactic potential fit function",
quad, ". Wait..."];
DiffFitGalPotPoints =
Table[Subscript[\[Phi], tot][GalPotPoints[[nr, 1]],
GalPotPoints[[nr, 2]]] - (quad /. X -> GalPotPoints[[nr, 1]] /=