Prev: Continued Fraction Trouble
Next: Misprint
From: Marco Masi on 18 Jul 2010 01:04
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,
" %"];
From: janos on 19 Jul 2010 02:09
On j=FAl. 18, 07:04, Marco Masi <marco.m...(a)ymail.com> 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 successful: however, I tried several functions (polynomial, logarithmic, 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 parameterizing 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 = 10000=
0;
>
> 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; b=
3 =
> = 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,
> " %"];

Something is wrong here:
DiffFitGalPotPoints =
Table[Subscript[\[Phi], tot][GalPotPoints[[nr, 1]],
GalPotPoints[[nr, 2]]] - (quad /. X -> GalPotPoints[[nr, 1]] /=
..
Y -> GalPotPoints[[nr, 2]]), {nr, 1, Length[GalPotPoints]}]=
;
Would you correct it please.

J=E1nos

From: Kevin J. McCann on 20 Jul 2010 03:45
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,
> " %"];

 | 
Pages: 1
Prev: Continued Fraction Trouble
Next: Misprint