Message Boards

WOLFRAM COMMUNITY

6692 Views

5 Replies

3 Total Likes

View groups...

Follow this post

Share this post:

GROUPS:

How can I use Manipulate with a NIntegrate defined function?

Eliseo Pavone

Posted 7 years ago

Attachments:

POSTED BY: Eliseo Pavone

5 Replies

Sort By:

Michael Rogers

Michael Rogers, Emory University

Posted 7 years ago

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[I2 Pi/nlines]# &, NestList[z0# &, 23., npts - 1], nlines - 1] ], "BoundaryElements" -> Join[ {LineElement[{{1, 2}, {2, 3}, {3, 4}, {4, 1}}]}, Table[LineElement[Partition[4 + knpts + 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]; 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[7Pi/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.01394106510^5interp[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.87241676710^3Exp[(-r)/2.676]Exp[(1)] * Integrate[Exp[-Abs[z]], {z, -Infinity, +Infinity}] /. r -> Last@array}, None} ] ]

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

Gianluca Gorni

Gianluca Gorni, University of Udine

Posted 7 years ago

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, (425^2)^(1/4), 1/100}]; interp = Interpolation[tbl, InterpolationOrder -> 1]; g[x_?NumericQ, y_?NumericQ, t_?NumericQ] = 2.87241676710^3Exp[(-Sqrt[x^2 + y^2])/2.676] Exp[(-Abs[z] + 1 - Cos[2(ArcTan[y/x] - t) + (2/(Tan[7Pi/90]))* Log[Sqrt[x^2 + y^2]/4.9]])]Exp[Abs[z]] Integrate[Exp[-Abs[z]], {z, -Infinity, +Infinity}] + 7.01394106510^5interp[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]

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

Michael Rogers

Michael Rogers, Emory University

Posted 7 years ago

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

Eliseo Pavone

Posted 7 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

Gianluca Gorni

Gianluca Gorni, University of Udine

Posted 7 years ago

Have you tried `Compile`?

POSTED BY: Gianluca Gorni

Reply to this discussion

Reply Preview

Attachments

Remove Add a file to this post

Follow this discussion

or Discard

Group Abstract

Feedback