I am a beginner using Mathematica. I am doing a physics experiment that involves this 3D surface; the function $ f(x, y) $. I need a program that can give me the contour level in the $ z $-axis which encloses an area that I specify.

This is a 3d view of the function

So far I have managed to create an algorithm that can accept a contour level and then output the area that it encloses:

The blue region is the area that is calculated

Here is the code: (credits to Kuba :Area between Contours in ContourPlot)

f[x_, y_] := 
 1/0.005*2.86*10^(-5)*
  Sin[Pi/2 - 
    ArcTan[(0.04367/35)*((Sin[2*Pi*35*y])^2)*2*(Pi/0.00763)*
      Cos[2*(Pi/0.00763)*x]]]*(2*
     Pi*35*(0.04367/35)*((Sin[2*Pi*35*y])^2)*Sin[2*(Pi/0.00763)*x]*
     Cos[2*Pi*35*
       y] + (2*9.81*(0.023 - 
         2*1.9*10^(-3) - (0.04367/35)*((Sin[2*Pi*35*y])^2)*
          Sin[2*Pi/0.00763*x]))^0.5)

plot = RegionPlot[0.0036 <= f[x, y], {x, 0, 0.007632}, {y, 0, 1/35}, 
  PlotPoints -> 100]

poly = Cases[Normal@plot, Polygon[n_] :> n, \[Infinity]]
Graphics`Mesh`MeshInit[];
PolygonArea /@ poly // Total

The value 0.0036, seen in RegionPlot, is the contour level in this case. And the program outputs an area value of 0.0000394488 when I run it. But how do I specify an area value and then get the contour level out of the program?

This is for my school project which is due soon; any help is immensely appreciated!!!

New contributor
Steve Gordan is a new contributor to this site. Take care in asking for clarification, commenting, and answering. Check out our Code of Conduct.
  • The solution in the question was originally written by Kuba here. – C. E. 2 days ago
  • Hi, thanks for the reply! But I referred to Kuba's code in writing my code above, but its essentially the same thing; they go from a specified contour level to the area, but for my experiment, I need to find the contour level that has a specified area enclosed in it. So its the reverse process. – Steve Gordan 2 days ago
  • I am so sorry Sir! I had no idea. – Steve Gordan 2 days ago
  • Write a function area[level] that returns the area of f[x, y] > level and use FindRoot[] to find the desired level – Michael E2 2 days ago
  • @MichaelE2 Thank you! However, I am terrible at coding. Please guide me further. This is what I have so far : f[x_, y_] := ..... area[level_] = RegionPlot[level <= f[x, y], {x, 0, 0.007632}, {y, 0, 1/35}, PlotPoints -> 100] poly = Cases[Normal@plot, Polygon[n_] :> n, [Infinity]] GraphicsMeshMeshInit[]; PolygonArea /@ poly // Total FindRoot[area[0.000039], {x = 0}] – Steve Gordan yesterday

As suggested by Micheal E2, use FindRoot

$Version

(* "11.3.0 for Mac OS X x86 (64-bit) (March 7, 2018)" *)

Clear["Global`*"]

f[x_, y_] := 
 Evaluate@Simplify[
   Rationalize[
    1/0.005*2.86*10^(-5)*
     Sin[Pi/2 - 
       ArcTan[(0.04367/35)*((Sin[2*Pi*35*y])^2)*2*(Pi/0.00763)*
         Cos[2*(Pi/0.00763)*x]]]*(2*Pi*35*(0.04367/35)*((Sin[2*Pi*35*y])^2)*
        Sin[2*(Pi/0.00763)*x]*
        Cos[2*Pi*35*
          y] + (2*9.81*(0.023 - 
            2*1.9*10^(-3) - (0.04367/35)*((Sin[2*Pi*35*y])^2)*
             Sin[2*Pi/0.00763*x]))^0.5), 0]]

{xmin, xmax} = {0, 477/62500};
{ymin, ymax} = {0, 1/35};
{fmin, fmax} = #[{f[x, y], xmin <= x <= xmax, ymin <= y <= ymax}, {x, 
     y}] & /@ {NMinValue, NMaxValue}

(* {0.00244846, 0.00419155} *)

area[level_?NumericQ] := Module[{plot, polys},
  plot = RegionPlot[level <= f[x, y], {x, xmin, xmax}, {y, ymin, ymax}, 
    PlotPoints -> 100];
  polys = Cases[Normal@plot, _Polygon, ∞];
  RegionMeasure /@ polys // Total]

Checking with the level from your example

area1 = area[0.0036]

(* 0.0000394488 *)

Using FindRoot to reverse the process

threshold[rgnMeasure_?(0 < # < (xmax - xmin) (ymax - ymin) &)] := 
 FindRoot[area[level] == rgnMeasure, {level, (fmax + fmin)/2}]

Unfortunately, this is quite slow

threshold[area1] // AbsoluteTiming

(* {110.345, {level -> 0.0036}} *)
  • OHH YESSSS THANK YOU SOO MUCHHHHHHHHHH!!!!!!!!!!!! I am forever in your debt, Mr Bob Hanlon. – Steve Gordan yesterday
  • Hi! Is there a way to decrease the precision or accuracy of the calculation to speed it up? Because I only need 4 significant figures – Steve Gordan 11 hours ago
  • @SteveGordan - In the RegionPlot replace the option PlotPoints -> 100 with MaxRecursion -> 3. That should cut the time to about 25% of the original time. If that is not enough, play with the options to FindRoot – Bob Hanlon 10 hours ago
  • Thanks! It's much faster now :))) – Steve Gordan 9 hours ago

Let me just add my five cents on this problem. The idea that I pursued was the same as Bob but I tried to make it faster by using one of the Monte Carlo methods of NIntegrate instead of RegionPlot.

This is my function for finding the area in terms of level threshold:

area[thresh_] := NIntegrate[
  Boole[f[x, y] >= thresh],
  {x, 0, xmax},
  {y, 0, ymax},
  Method -> "AdaptiveMonteCarlo"
  ]

area[0.0036] // AbsoluteTiming

{0.030328, 0.0000390467}

This is quite fast and the result is almost the same as the one given by OP's method. Unfortunately, I found that FindRoot, NMinimize etc. do not play nice with this function, I'm guessing because it is stochastic. The Jacobian is singular, FindRoot says, for example. Using a similarly stochastic method of NMinimize such as SimulatedAnnealing also didn't work, it is still unclear to me why.

So I had to resort to a workaround which rather put a blemish on the whole thing, here it is:

table = Table[{Quiet@area[x], x}, {x, 0.003, 0.004, 0.001 0.003}];
areaInTermsOfLevel = Interpolation[table];
areaInTermsOfLevel[0.00003930850798691388`]

0.00360097

This works but yeah, it's not very pretty, and also evaluating area a thousand times is pretty inefficient. Another problem is finding the limits 0.003 and 0.004. The reason is that e.g. the largest level allowed is too large (0.005 say) then the area will be 0 for multiple values, and that will cause problems for Interpolation. (Bob's answer shows how to get the thresholds that should work, in theory. The smallest and largest values of f[x, y])

  • THANK YOUU!!! What do I have to do achieve your coding prowess? Please tell me. The way you throw this jargon around is mesmerizing. I want to be like you. What's your secret? What did you study? – Steve Gordan yesterday
  • @SteveGordan I mostly learned by reading and contributing here on Stack Exchange, a well as by studying physics at the university. – C. E. yesterday

Your Answer

Steve Gordan is a new contributor. Be nice, and check out our Code of Conduct.

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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