Group Abstract Group Abstract

Message Boards Message Boards

How can I use Manipulate with a NIntegrate defined function?

Eliseo Pavone
Posted 8 years ago

Hi, I'm trying to model a spiral galaxy using some density distributions and display them with a 2D density plot. The functions are integrated along a line of sight in order to reduce the 3 variables functions to 2 variables function. In order to make a rotating plot I'm trying to use Manipulate. Here it is my code, which gives me as output only a static image that doesn't rotate moving the slider. The code is this:

LogarithmicScaling[x_, min_, max_] := Log[x/min]/Log[max/min]
f[x_?NumericQ, y_?NumericQ, t_?NumericQ] := 
 NIntegrate[
  2.872416767*10^3*Exp[(-Sqrt[x^2 + y^2])/2.676]*Exp[(-Abs[z] + 1- Cos[2*(ArcTan[y/x] - t) + (2/(Tan[7*Pi/90]))*Log[Sqrt[x^2 + y^2]/4.9]])] +7.013941065*10^5*Exp[-Sqrt[x^2 + y^2 + z^2]/0.09866] + (7.321354473)/(1 + (Sqrt[x^2 + y^2 + z^2]/3.936)^2)^(5/2), {z, -Infinity, +Infinity}]

plotter[min_, max_, NumberOfTicks_] := 
 Manipulate[DensityPlot[f[x, y, t], {x, -25, 25}, {y, -25, 25}, 
   PlotRange -> Full, ColorFunctionScaling -> False, 
   ColorFunction -> (ColorData["SunsetColors"][
       LogarithmicScaling[#, min, max]] &), PlotPoints -> 50, 
   PlotLegends -> 
    BarLegend[{ColorData["SunsetColors"], {0, 1}}, 
     LegendMarkerSize -> 370, 
     Ticks -> ({LogarithmicScaling[#, min, max], 
          ScientificForm[#, 1]} & /@ (min (max/min)^
           Range[0, 1, 1/NumberOfTicks]))]], {t, 0, 100}, 
  AnimationRunning -> False, SynchronousUpdating -> False ]
plotter[10^1, 10^5, 4]

The static image which has to rotate is shown in the attachment. If I plot the function fixing t with different values, it actually rotates, but the manipulate function seems to not work at all. Any idea?

Attachments:
POSTED BY: Eliseo Pavone
5 Replies
Sort By:

Building on Gianluca's ideas, one can do two things, separate out further the smallest bit of the computation that depends on the dynamic variable t and manually construct a mesh that efficiently represents the plot (rather than use DensityPlot's rectangular mesh, which is rather awkward for this plot).

Needs["NDSolve`FEM`"];

rate = 3;  (* can set it to 2 for a smaller mesh *)
npts = rate * 40; (* number of steps in the spiral lines *)
nlines = 15;  (* can set it to 10 for a smaller mesh *)
z0 = Exp[(-0.1 + 0.5 I )/rate];
bmesh = ToBoundaryMesh["Coordinates" -> Join[
     {{-25., -25.}, {25., -25.}, {25., 25.}, {-25., 25.}},
     ReIm@
      Flatten@NestList[Exp[I*2 Pi/nlines]*# &, 
        NestList[z0*# &, 23., npts - 1], nlines - 1]
     ],
   "BoundaryElements" -> Join[
     {LineElement[{{1, 2}, {2, 3}, {3, 4}, {4, 1}}]}, 
     Table[LineElement[Partition[4 + k*npts + Range@npts, 2, 1]],
      {k, 0, nlines - 1}]
     ]
   ];
emesh = ToElementMesh[bmesh]
Show[
 emesh["Wireframe"],
 emesh["Wireframe"["MeshElement" -> "BoundaryElements",  "MeshElementStyle" -> Red]]
 ]
gc = First@Cases[emesh["Wireframe"], _GraphicsComplex, Infinity];

enter image description here

It has almost 8000 points. You probably don't want it to get too much above 10000.

emesh["Coordinates"] // Length
(*  7817  *)

Then

Clear[g0];
g0[t_] := Exp[(-Cos[2*(a - t) + (2/(Tan[7*Pi/90]))*Log[r/4.9]])];

With[{array = {ArcTan[#[[1]], #[[2]]], Sqrt[#[[1]]^2 + #[[2]]^2]} &@
    Transpose@First@gc},
 Manipulate[
  Graphics[
   Append[
    gc,
    VertexColors -> 
     Dynamic[ColorData["SunsetColors"] /@ 
       LogarithmicScaling[(scale*(g0[t] /. Thread[{a, r} -> array]) + 
          offset), 10.^1, 10.^5
        ]]
    ]
   ],
  {t, 0, 10},
  {{offset, +7.013941065*10^5*interp[r^2] + 
      Integrate[(7.321354473)/(1 + (Sqrt[r^2 + z^2]/3.936)^2)^(5/2), {z, -Infinity, +Infinity},
        Assumptions -> r^2 >= 0] /. 
     Thread[{a, r} -> array]}, None},
  {{scale, 
    2.872416767*10^3*Exp[(-r)/2.676]*Exp[(1)] * Integrate[Exp[-Abs[z]], {z, -Infinity, +Infinity}] /.
       r -> Last@array}, None}
  ]
 ]

enter image description here
POSTED BY: Michael Rogers

Your integration is too onerous. Two of the summands have simple closed forms, while the third can be pre-computed and interpolated for speed. This is a version that is responsive enough:

tbl = Table[{u^4, 
    NIntegrate[
     Exp[-Sqrt[u^4 + z^2]/0.09866], {z, -Infinity, +Infinity}]},
   {u, 0, (4*25^2)^(1/4), 1/100}];
interp = Interpolation[tbl, InterpolationOrder -> 1];
g[x_?NumericQ, y_?NumericQ, t_?NumericQ] = 
  2.872416767*10^3*Exp[(-Sqrt[x^2 + y^2])/2.676]*
    Exp[(-Abs[z] + 1 - 
       Cos[2*(ArcTan[y/x] - t) + (2/(Tan[7*Pi/90]))*
          Log[Sqrt[x^2 + y^2]/4.9]])]*Exp[Abs[z]]*
    Integrate[Exp[-Abs[z]], {z, -Infinity, +Infinity}] +
   7.013941065*10^5*interp[x^2 + y^2] + 
   Integrate[(7.321354473)/(1 + (Sqrt[x^2 + y^2 + z^2]/3.936)^2)^(5/
        2), {z, -Infinity, +Infinity}, 
    Assumptions -> x^2 + y^2 >= 0];
plotter[min_, max_, NumberOfTicks_] := 
  Manipulate[
   DensityPlot[g[x, y, t], {x, -25, 25}, {y, -25, 25}, 
    PlotPoints -> 50],
   {t, 0, 100}];
plotter[10^1, 10^5, 4]
POSTED BY: Gianluca Gorni

Have you tried Compile?
POSTED BY: Gianluca Gorni
Eliseo Pavone
Eliseo Pavone
Posted 8 years ago

Thank you Gianluca, your advice has been really helpful. At the end I chose to do a table with all the plots and create an .avi file because it was still too slow with manipulate.
POSTED BY: Eliseo Pavone

For me it works if I fix the syntax error with BarLegend not having the ] in the right place. However, it takes a long time to compute. I assume it rotated eventually when I moved t, but I wasn't looking when the computation finished.
POSTED BY: Michael Rogers
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard