Finding area under a region of a Parametric Plot using mesh functions

I have the following Mathematica code:

\$Assumptions = {r \[Element] Reals, r >= 0, rh \[Element] Reals,
rh > 0, g \[Element] Reals, g > 0, l \[Element] Reals, l > 0};
M[rh_, g_] := (1 + rh^2/l^2) (rh^2 + g^2)^(3/2)/(2 rh^2)
f[r_, rh_, g_] := 1 - (2 M[rh, g] r^2)/(g^2 + r^2)^(3/2) + r^2/l^2
T[rh_, g_] := D[f[r, rh, g], r]/(4 \[Pi]) /. r -> rh
Veff[r_] := f[r, rh, g]*(L^2/r^2)
r0 = Solve[Veff'[r] == 0, r, Reals][[2]]
\[Lambda][rh_,
g_] = (Sqrt[-((r^2 f[r, rh, g])/L^2) Veff''[r]]) /. r0 //
Simplify;
l = 1;
(tp1 = {#[[1]], \[Lambda][rh1 = (rh /. #[[2]]), 0.02]} &@
Maximize[{T[rh, 0.02], 0.01 < rh < 1}, rh] // FullSimplify) // N;
(tp2 = {#[[1]], \[Lambda][rh2 = (rh /. #[[2]]), 0.02]} &@
Minimize[{T[rh, 0.02], rh1 < rh < 1}, rh] // FullSimplify) // N

pplt = ParametricPlot[{T[rh, 0.02], \[Lambda][rh, 0.02]}, {rh, 0.01,
8}, PlotRange -> {{0, 1}, {0, 10}}, AspectRatio -> 0.5,
PlotPoints -> 50, ImageSize -> 800,
PlotTheme -> {"Scientific", "BoldColor"},
FrameLabel -> {"Temperature \[Rule]", "Lyapunov Exponent \[Rule]"},
LabelStyle -> {Bold, Black, 22}, PlotStyle -> Thickness[0.004],
ColorFunction ->
Function[{T, \[Lambda], rh},
If[rh <= rh1, Blue, If[rh <= rh2, Red, Green]]],
ColorFunctionScaling -> False,

Prolog -> {{Gray, Dashing[0.008], Line[{tp1, {tp1[[1]], 0}}],
Line[{tp2, {tp2[[1]], 0}}]}, {Darker[Black],
Line[{{((tp1 + tp2)/2)[[1]], 0.5}, {((tp1 + tp2)/2)[[1]],
9}}]}, {Gray, Thick, Dashing[0.001], Circle[tp1, {0.01, .2}],
Circle[tp2, {0.01, .2}]}}];

P3 = Legended[pplt,
Placed[SwatchLegend[{Blue, Red, Green}, {"Small BH",
"Intermediate BH", "Large BH"}, LabelStyle -> {Bold, Black, 20},
LegendMarkers -> "SphereBubble", LegendMarkerSize -> 20,
LegendFunction -> (Framed[#, RoundingRadius -> 10,
FrameStyle -> None, Background -> GrayLevel[0.95]] &),
LegendLabel ->
Placed["Black Holes", Left,
Rotate[Style[#, 16], 90 Degree] &]], {0.2, .6}]]


which gives me the following result:

How can I find the area under Region 1 and Region 2 using Mesh functions ? I also wish to find the position of the line (Basically the value of T) for which the Region 1 and Region 2 have equal areas. What approach should I use for that, without succumbing to trial-and-error methods?

• What does the restriction "using mesh functions" mean? Commented Aug 25, 2023 at 7:08
• We can also write λ as a function of T, by substitution which can be used to find the numerical integration using proper limits. This is quite a cumbersome process. I was looking for a way, where I can find the specified area, directly from the plot itself. I have previously seen questions where they talk about finding area using mesh regions in a different context like -> mathematica.stackexchange.com/a/119220/78049, So I thought of using such approach for my problem as well. Commented Aug 25, 2023 at 7:32
• @codepr Perhaps my answer is what you're looking for? Commented Aug 25, 2023 at 7:36
• @UlrichNeumann yes it works. Thanks a lot for this! Actually, my main objective was to find the black line which intersects the whole curve into two equal areas, or in other words the value of T, for which a perpendicular line when drawn through it justifies Maxwell construction(en.wikipedia.org/wiki/Maxwell_construction), I think this requires a different question all together. Again, really appreciate your help ! Commented Aug 25, 2023 at 7:42
• See my modified answer Commented Aug 25, 2023 at 8:16

Hope I understood your question correctly.

Using P3 as start, first take all the points

pts = Cases[Normal[P3], Line[p_] :> p, -1][[1]];


Now take three subsets defined by range tp1,tp2

blue = Select[pts, #[[1]] <= tp1[[1]] && #[[2]] >= tp1[[2]] &];
red = Select[pts,tp2[[1]] <= #[[1]] <= tp1[[1]] && tp2[[2]] <= #[[2]] <= tp1[[2]] &] ;
green = Select[pts ,  #[[1]] >=   tp2[[1]] &&   #[[2]] <= tp2[[2]] &] ;


Interpolate these points

ipblue = Interpolation[blue];
ipred = Interpolation[red];
ipgreen = Interpolation[green];


and integrate

{Ablue, Ared, Agreen} =NIntegrate[{ipblue[tp], ipred[tp], ipgreen[tp]}, {tp,tp2[[1]],tp1[[1]]}] // Quiet


The two areas you 're looking for follow to

{Ablue - Ared, Ared - Agreen}
(*{2.82274, 1.9862}*)


Hope it helps!

calculation of midpoint T

reg1[T_?NumericQ] :=NIntegrate[ ipblue[tp] - ipred[tp] , {tp, T, tp1[[1]]} ]
reg2[T_?NumericQ] := NIntegrate[  ipred[tp] - ipgreen[tp] , {tp, tp2[[1]], T} ]

line1 = Line[Map[{#, reg1[#]} &, Subdivide[tp2[[1]], tp1[[1]], 20]]];
line2 = Line[Map[{#, reg2[#]} &, Subdivide[tp2[[1]], tp1[[1]], 20]]];
pT=RegionIntersection[line1, line2](*Point[{{0.669968, 0.707078}}]*)

T=pT[[1, 1]][[1]] (*0.669968*)


More direct solution using FindRoot (NDSolve and NMinimize won't evaluate)

T =.
FindRoot[reg1[T] == reg2[T], {T, tp2[[1]], tp1[[1]]}]
(*{T -> 0.669965} *)

• Thanks a lot for the much needed help! Commented Aug 25, 2023 at 13:06
• @codebpr You're welcome Commented Aug 25, 2023 at 13:19

After answering your question, here I'll add a supplement which determines tp1,tp2 using MeshFunctions

plot = ParametricPlot[{T[rh, 0.02], \[Lambda][rh, 0.02]}, {rh, 0.01,
8}, PlotRange -> {{0, 1}, {0, 10}}, AspectRatio -> 0.5,
PlotPoints -> 50, ImageSize -> 800,
MeshFunctions ->
Function[{x, y, rh},
Evaluate[Derivative[1, 0][ T][rh, 1/50]/1(*Derivative[ 1,
0][\[Lambda]][rh,1/50]*)]], Mesh -> {{0}}]


get points

 pts = plot[[1]][[1, 1]]; (* all points*)
index = Cases[plot, Point[ p_] :>  p, -1][[1]]  (* index meshpoints*)
tp=pts[[index]] (*{tp1,tp2}*)
(*{{0.274999, 1.5824}, {0.951599, 7.78498}}*)

Show[{ListPlot[pts , PlotRange -> All],Graphics[{Red,PointSize[Large], Point[tp]}]}, PlotRange -> {{0, 1}, {0, 10}}]


• Thanks a lot for this :) Just one small doubt. Suppose I fix T=0.35 (which is basically related to transition temperature for this system) from my side and now wish to find the ratio of the new areas, what tweaks should I make in the previous code snippet of yours? Commented Aug 25, 2023 at 14:40
• Something like reg1[.35]/reg2[.35] Commented Aug 25, 2023 at 15:18