Re: Fitting complicate functions
- To: mathgroup at smc.vnet.net
- Subject: [mg111106] Re: Fitting complicate functions
- From: "Kevin J. McCann" <Kevin.McCann at umbc.edu>
- Date: Tue, 20 Jul 2010 03:45:09 -0400 (EDT)
- References: <i1u21r$7sd$1@smc.vnet.net>
It is possible that a sparse Bayesian curve fit (regression) may do it.
The idea of such a fit is to start with a large number of fitting
functions; however, it is found that most end up being zero. If you are
interested, I could look into it for your case.
Kevin
Marco Masi wrote:
> 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]] /=
> .
> Y -> GalPotPoints[[nr, 2]]), {nr, 1, Length[GalPotPoints]}]=
> ;
> MaxDiffPot =
> Max[Table[DiffFitGalPotPoints[[nr]], {nr, 1, Length[GalPotPoints]}]];
> MinDiffPot =
> Min[Table[DiffFitGalPotPoints[[nr]], {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]]
> Print["The above figure: plot of the difference between the Galactic \
> potential and the Galactic potential fit function."];
> Print["Minimum Galactic potential of ", MinPot,
> " Maximun positive difference= ", MaxDiffPot,
> " Minimum negative difference= ", MinDiffPot];
> PotErrRange = 100*(MaxDiffPot + Abs[MinDiffPot])/Abs[MinPot];
> Print["Maximum error of fit over Galactic potential= ", PotErrRange,
> " %"];