9
$\begingroup$

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}]

enter image description here

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.90873735713213132056690405656577392873885897791564608027113625396018946387876505698136804341595759`100 with error 0.0007620609074285774952076996909671179060090950094015269438194254831966603709711198623965417774800435432`100, 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?

Improve this question
$\endgroup$
2
  • 2
    $\begingroup$ Meshing is still machine-precision only, I believe. So you wouldn't succeed with mathods based on it. $\endgroup$ Commented Aug 9, 2024 at 14:24
  • 2
    $\begingroup$ 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 $\endgroup$ Commented Aug 9, 2024 at 17:56

2 Answers 2

Reset to default
11
$\begingroup$

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.089257904195311118302821843453118688905478368275606690887764085319790843430949003608616784986395099690694426729361191802957181759530004825591584159667415684061336910338923062197881745123171307854734395769307124116366790293855586961037970607233501172490718249413687613573868432825911159923776112844980253023512869693814921828936`300.; 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 *)
Improve this answer
$\endgroup$
10
$\begingroup$

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.84554907924043310148869075386010574213`20.412322424771425}], 0}, WorkingPrecision -> 200, MaxRecursion -> 1000]

The two results and their difference:

{2.08925790419531111830282184345311868890547836827560669088776408531979084343094900360861678498639509969069442672936119180295718175953000482559158415966741568406133691033892306219788174512317131234402512557826284126`200., 2.089257904195311118302821843453118688905478368275606690887764085319790843430949003608616784986395099690694426729361191802957181759530004825591584159667415684061336910338923062197881745123171307854734395769307124116366790293855586961037970607233501172490718249413687613573868432825911159923776112844980253023512869693814921828936`300. } // 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}]}} *)
Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.