Message Boards Message Boards

GROUPS:

How can I use Manipulate with a NIntegrate defined function?

Posted 11 months ago
1207 Views
|
5 Replies
|
3 Total Likes
|

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?

Attachment

Attachments:
5 Replies

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.

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 11 months 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.

Have you tried Compile?

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

Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract