# High-precision solution for the area of a region

I want to solve the area of the region ImplicitRegion[Exp[-(y + x^3)^2] > 1/32 y^2 + Exp[Sin[y]], {x, y}] with high precision (like 100 digits). The region looks like

ContourPlot[Exp[-(y + x^3)^2] == 1/32 y^2 + Exp[Sin[y]], {x, -2, 2}, {y, -4, 0}]


However, I've tried different methods to calculate the area, but these methods gave different values.

The first one is provided by Area

Area[ImplicitRegion[
Exp[-(y + x^3)^2] > 1/32 y^2 + Exp[Sin[y]], {x, y}],
WorkingPrecision -> 50, PrecisionGoal -> 50]
(*2.1150271808931174*)


We can see that PrecisionGoal and WorkingPrecision don't work here.

The second one is provided by NIntegrate using Boole

NIntegrate[
Boole[Exp[-(y + x^3)^2] > 1/32 y^2 + Exp[Sin[y]]], {x, -2,
2}, {y, -4, 0}, WorkingPrecision -> 50, PrecisionGoal -> 50]


It gives a lot of warning and returns 1.90873735713213132056690405656577392873885897791564608027113625396018946387876505698136804341595759100 with error 0.0007620609074285774952076996909671179060090950094015269438194254831966603709711198623965417774800435432100, which are absolutely wrong.

The third one is provided by NIntegrate using Element

NIntegrate[1, {x, y} \[Element]
ImplicitRegion[
Exp[-(y + x^3)^2] > 1/32 y^2 + Exp[Sin[y]], {{x, -2, 2}, {y, -4, 0}}],
WorkingPrecision -> 50, PrecisionGoal -> 50]


It returns 2.08924, and we can see that PrecisionGoal and WorkingPrecision don't work here.

The last one is using DiscretizeRegion before using Area

Area[DiscretizeRegion[
ImplicitRegion[
Exp[-(y + x^3)^2] >
1/32 y^2 + Exp[Sin[y]], {{x, -2, 2}, {y, -4, 0}}],
AccuracyGoal -> 50, MaxCellMeasure -> #], AccuracyGoal -> 50,
WorkingPrecision -> 50] & /@ {0.01, 0.0001, 0.00001}
(*{2.0892495973037697, 2.0892577898447136, 2.089257792018752}*)


Again, WorkingPrecision and AccuracyGoal don't work here, and MaxCellMeasure will affect the final result.

So my question is: how to obtain a high-precision result with mathematica?

• Meshing is still machine-precision only, I believe. So you wouldn't succeed with mathods based on it. Commented Aug 9 at 14:24
• The region shown by your ContourPlot is not the only one were the inequality holds. Exp[-(y + x^3)^2] > 1/32 y^2 + Exp[Sin[y]] /. {x->-1.6,y->4.1} gives True Commented Aug 9 at 17:56

Here's another way to find the area of the loop below the x-axis, using the Green's-theorem area formula and NDSolve to perform the line integral. The value for area300 is the 300-digit result from my other answer, used to check the accuracy of the NDSolve result. For some reason it takes just a few seconds to get a 100-digit-accurate result, but it takes 100 sec. with WorkingPrecision -> 300, PrecisionGoal -> 125 to yield 127 accurate digits.

area300 =  2.089257904195311118302821843453118688905478368275606690887764085319790843430949003608616784986395099690694426729361191802957181759530004825591584159667415684061336910338923062197881745123171307854734395769307124116366790293855586961037970607233501172490718249413687613573868432825911159923776112844980253023512869693814921828936300.;
eqn = Exp[-(y + x^3)^2] == 1/32 y^2 + Exp[Sin[y]];
odeTraj =
D[{x[t], y[t]},
t] == -Cross[
D[eqn /. Equal -> Subtract, {{x, y}}] /. {x -> x[t],
y -> y[t]}];
icTraj = {x[0], y[0]} == {0, 0};
odeArea = A'[t] == 1/2 Det[NestList[D[#, t] &, {x[t], y[t]}, 1]];
icArea = A[0] == 0;
NDSolveValue[{odeTraj, odeArea, icTraj, icArea,
WhenEvent[x[t] == 0 && Abs[y[t]] < 1/100, "StopIntegration"]}
, A[Indexed[A@"Domain", {1, 2}]]
, {t, 0, Infinity}
, Method -> "Extrapolation"
, WorkingPrecision -> 250, PrecisionGoal -> 100] // AbsoluteTiming
Last@% - area300 // N
(*
{6.33992,
2.089257904195311118302821843453118688905478368275606690887764085319790843430949003608616784986395099552233052235378808949201217589524587321214765400244628245451592405564091800205665874730607674476527006975559216817071559694893332831170747468620947586}

-1.38461*10^-100
*)


You can solve the equation fo $$x$$ fairly easily. Then you can NIntegrate at working precisions 200 and 300, which have precision goals of 100 and 150 respectively. Their difference gives an estimate of the error of the first integral. It takes a minute or so.

NIntegrate[-Surd[-y - Sqrt[Log[32/(32 E^Sin[y] + y^2)]], 3] +
Surd[-y + Sqrt[Log[32/(32 E^Sin[y] + y^2)]], 3],
{y,
Root[{-32 + 32 E^Sin[#1] + #1^2 &,
-2.8455490792404331014886907538601057421320.412322424771425}],
0},
WorkingPrecision -> 200, MaxRecursion -> 1000]


The two results and their difference:

{2.08925790419531111830282184345311868890547836827560669088776408531979084343094900360861678498639509969069442672936119180295718175953000482559158415966741568406133691033892306219788174512317131234402512557826284126200.,
2.089257904195311118302821843453118688905478368275606690887764085319790843430949003608616784986395099690694426729361191802957181759530004825591584159667415684061336910338923062197881745123171307854734395769307124116366790293855586961037970607233501172490718249413687613573868432825911159923776112844980253023512869693814921828936300.
} // Differences

(*  {-4.4892907*10^-192}  *)


You can solve for the interval integration thus:

Solve[Log[32/(32 E^Sin[y] + y^2)] == 0 && -3 < y <= 0, {y}]
(*
{{y -> 0}, {y ->
Root[{-32 + 32 E^Sin[#1] + #1^2 &, -2.8455490792404331015}]}}
*)