2
$\begingroup$

I have to use ContourPlot with a complicated function depending on 2 parametrs (that I cannot report here) that contains numerical integrations.

Here is a simple example to clarify my problem:

f[a_?NumberQ, b_?NumberQ] := NIntegrate[a*x + b, {x, 0, 1}];
ContourPlot[{f[a, b] == 3, f[a, b] == 6, f[a, b] == 10}, {a, 0, 
   10}, {b, 0, 10}] // Timing

enter image description here

As you can see, since Nintegrate is evaluated in every point, the computation is very slow (over 3 seconds on my laptop, even for this simple case): how can I optimize this?

$\endgroup$
  • 2
    $\begingroup$ ContourPlot[f[a, b], {a, 0, 10}, {b, 0, 10}, Contours -> {3, 6, 10}] // Timing isn't probably the best solution, but it shaves the time by 40% $\endgroup$ – Dr. belisarius Feb 10 '15 at 17:00
  • 1
    $\begingroup$ Maybe also play with the number of PlotPoints. $\endgroup$ – Jens Feb 10 '15 at 17:34
  • 2
    $\begingroup$ this is redundantly computing the Integral at some points, saving values ( f[a_?NumberQ, b_?NumberQ] := f[a,b]= NIntegrate.. ) cuts the time in half for this example. $\endgroup$ – george2079 Feb 10 '15 at 20:51
5
$\begingroup$

Update: Using the memoization trick suggested by @george2079 in the comments combined with MeshFunctions we get the same picture in 0.015625 seconds:

fa[a_?NumericQ, b_?NumericQ] := fa[a, b] = NIntegrate[a*x + b, {x, 0, 1}];
First@Timing[ContourPlot[fa[a, b], {a, 0, 10}, {b, 0, 10}, Contours -> {}, 
    ContourShading -> None, MeshFunctions -> {#3 &}, 
    Mesh -> {{{3, Red}, {6, Green}, {10, Blue}}}, 
    MeshStyle -> Thick] /. Polygon[__] -> Sequence[]]
(* 0.015625 *)

Original post:

Using MeshFunctions with no contours

First@Timing[cp1 = ContourPlot[f[a, b], {a, 0, 10}, {b, 0, 10}, 
   Contours -> {}, ContourShading -> None,
   MeshFunctions -> {#3 &}, Mesh -> {{{3, Red}, {6, Green}, {10, Blue}}}, 
   MeshStyle -> Thick]/. Polygon[__] -> Sequence[]]
(* 0.578125 *)

makes it faster

cp1

enter image description here

versus

First@Timing[ContourPlot[f[a, b], {a, 0, 10}, {b, 0, 10}, Contours -> {3, 6, 10}]]
(* 1.625000 *)
$\endgroup$
  • $\begingroup$ I cannot avoid numerical integrations $\endgroup$ – mattiav27 Feb 10 '15 at 17:18
  • $\begingroup$ @belisarius, mattiav27, sorry i missed the key requirement of the question:) $\endgroup$ – kglr Feb 10 '15 at 17:20
  • $\begingroup$ Much better now. +1 cp1 /. GrayLevel[_] :> GrayLevel[1] $\endgroup$ – Dr. belisarius Feb 10 '15 at 17:46
  • $\begingroup$ It is important to see why this is fast: here the function value is computed on a regular grid and the contours are based on simple interpolation. In the approach in the question ContourPlot uses a root finding algorithm to accurately find the locus of the contours. $\endgroup$ – george2079 Feb 10 '15 at 21:19
  • $\begingroup$ @george, good point, thank you. Actually, MeshFunctions combined with your memoization suggestion cuts the timing to almost zero. $\endgroup$ – kglr Feb 10 '15 at 21:35
1
$\begingroup$

Here is another alternative for how to speed up the plotting if you want to retain all the functionality and options of ContourPlot: use ParallelTable:

f[a_?NumericQ, b_?NumericQ] := NIntegrate[a*x + b, {x, 0, 1}];

Timing[Show[
  ParallelTable[
   ContourPlot[f[a, b], {a, 0, 10}, {b, 0, 10}, Contours -> {i}, 
    ContourShading -> None, 
    ContourStyle -> Hue[i/10]], {i, {3, 6, 10}}]]]

This is a little slower than kguler's solution, but I think it will depend on the actual function to be plotted, and on the complexity of the contour lines you want to see. My solution simply offloads the computation of each individual contour onto a separate kernel, by making a ParallelTable of ContourPlots for one contour each. Then I combine them with Show at the end.

You could in principle combine kguler's and my suggestion to speed things up still more, but I just want to point out the parallelization aspect in making the plot.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.