Evaluation order in DensityPlot

Question

I have a module that, when given a value of the variable a, returns a function of the variable x. My actual problem is longer and more complicated; here I merely provide an MWE, say:

TestModule[a_] := 
 Module[{}, {vals, sols} = 
   NDEigensystem[{D[\[Phi][u], {u, 2}] + 1/u D[\[Phi][u], {u, 1}] - 
      1/u^2 \[Phi][u] - 1/4 a^2 u^2 \[Phi][u], 
     DirichletCondition[\[Phi][u] == 0, u == 1]}, \[Phi][u], {u, 0, 
     1}, 5, Method -> {"PDEDiscretization" -> {"FiniteElement", 
        "MeshOptions" -> {"MaxCellMeasure" -> 10^-3}}}]; Return[\!\(
\*SubsuperscriptBox[\(\[Sum]\), \(i = 1\), \(5\)]\(
\*FractionBox[\(0.01\), \(
\*SuperscriptBox[\((x - vals[\([\)\(i\)\(]\)])\), \(2\)] + 
\*SuperscriptBox[\(0.01\), \(2\)]\)] sols[\([i]\)]\)\) /. {u -> 0}, 
   Module];]

For instance (just to make things concrete), TestModule[0.5] returns

2.51604*10^-6/(0.0001 + (14.7028 + x)^2) + 5.4659*10^-6/(
 0.0001 + (49.2393 + x)^2) + 8.92959*10^-6/(
 0.0001 + (103.52 + x)^2) - 0.0000128153/(
 0.0001 + (177.542 + x)^2) + 0.0000170653/(0.0001 + (271.302 + x)^2)

I would now like to make a DensityPlot in the variables a and x. However, if I write

DensityPlot[Evaluate@TestModule[a],{a,0,1},{x,-500,0}]

I get the error "The PDE coefficient 1/u^2+(a^2\u^2)/4 does not evaluate to a numeric scalar at the coordinate {0.5}; it evaluated to 4. +0.0625` a^2 instead". Naively it seems to me that DensityPlot is assuming I have an explicit expression in terms of a and then trying to evaluate the function, but this obviously fails at the level of DEigensystem. Instead, I would like that for each input of a that DensityPlot feeds the module, a function of x is obtained and then evaluated at the different values of x. Is there a way to ensure this and let DensityPlot first obtain the function of x for each given a, and only then evaluate at x? Please note that my actual problem is much more complicated and I would like to only evaluate the module once for each given a. I guess it all boils down to somehow modifying the evaluation order of DensityPlot. So ideally I'd like it to focus individually on each vertical slice, and then perform many of these slices corresponding to different values of a. Hopefully what I'm asking makes sense. Thank you in advance!

https://mathematica.stackexchange.com/a/26037 might be helpful here. — eyorble, Dec 10 '22 at 22:08
Thanks for your reply! I did not come across that post before and I'll read through everything else as well as it seems very useful. Writing TestModule[a_?NumericQ] in the definition does give a working DensityPlot. The problem is that then that the function gets evaluated at every pair (a,x). I really need it to only evaluated once at each a and then evaluate at x, otherwise the plot in my actual problem will take literally days to complete. Perhaps this is not possible with DensityPlot itself and I have to find an alternative. — Maisenberg, Dec 10 '22 at 22:57

score 0 · Accepted Answer · answered Dec 10 '22 at 23:44

The given system doesn't actually show much variation over the range given. Let us rewrite the system to allow memoization:

TestModule[a_] := TestModule[a] =
   Module[{}, {vals, sols} = 
       NDEigensystem[{D[φ[u], {u, 2}] + 1/u D[φ[u], {u, 1}] - 
             1/u^2 φ[u] - 1/4 a^2 u^2 φ[u], 
           DirichletCondition[φ[u] == 0, u == 1]}, φ[u], {u, 0, 1}, 5, 
     Method -> {"PDEDiscretization" -> {"FiniteElement", 
                 "MeshOptions" -> {"MaxCellMeasure" -> 10^-3}}}];
   Evaluate[Total[0.01/((# - vals)^2 + 0.01^2) (sols /. u -> 0)]] &]

Note that the return format is generalized slightly from the previous Sum notation, removed from the Return statement, and is returned as a function rather than an expression containing x. As an example:

TestModule[0.5]

2.51604*10^-6/(0.0001 + (14.7028 + #1)^2) + 5.4659*10^-6/( 0.0001 + (49.2393 + #1)^2) + 8.92959*10^-6/( 0.0001 + (103.52 + #1)^2) - 0.0000128153/( 0.0001 + (177.542 + #1)^2) + 0.0000170653/( 0.0001 + (271.302 + #1)^2) &

Specific x points can be extracted from this result as such: TestModule[0.5][x].

This allows the evaluation of a once for a large range of x. For example, after clearing Global`*, evaluation of Plot[TestModule[0.5][x], {x, -500, 0}] takes about a tenth of a second on my machine, resulting in:

.

It is also possible to generate a density plot by manually sampling the various a values, e.g.:

AbsoluteTiming[dat = Table[{a, x, TestModule[a][x]}, {a, 0, 1, 0.1}, {x, -500, 0, 0.5}];]

{0.616687, Null}

Relatively quick on my machine, for the number of points being generated here. Since the peaks are relatively narrow in x, the small step in x is chosen for ease of use. Nicer sampling techniques could probably use quite a smaller number of points though.

This can be plotted using ListDensityPlot:

ldp = ListDensityPlot[Flatten[dat, 1]];
ldp

There is not much variation in the result observed over the small range of a.

Thank you for this! Not much variation in a is perfectly possible since I just made up this precise differential equation for the MWE :) — Maisenberg, Dec 11 '22 at 13:57

Evaluation order in DensityPlot

1 Answers1