6
$\begingroup$

If I have the two following surfaces:

f[x_, y_] := y^3 - x^3 + 3 x^2 - 6 x - 4 y + 4
g[x_, y_] := x^3 + y^3 + x - 2 y - 1

And the two lines defined by the contours:

f[x,y] == 0
g[x,y] == 0

Then how can I get the region between the two contours?

I have the points of intersection which are the results in the first quadrant given by:

solution = {x, y} /. NSolve[ f[x, y] == g[x, y] == 0, Reals]
$\endgroup$
12
$\begingroup$

Here's an approach that involves some manual inspection, but works nonetheless.

First let's look at the two curves:

f[x_, y_] := y^3 - x^3 + 3 x^2 - 6 x - 4 y + 4
g[x_, y_] := x^3 + y^3 + x - 2 y - 1

ContourPlot[{f[x, y] == 0, g[x, y] == 0}, {x, -3, 3}, {y, -3, 3}]

enter image description here

We see there's 2 components we're interested in. To find them, let's look to see where f and g are negative:

RegionPlot[{f[x, y] < 0, g[x, y] < 0}, {x, -3, 3}, {y, -3, 3}]

enter image description here

So by inspection, we can see the region of interest is

f[x, y] < 0 && g[x, y] < 0 && y > 0 (* top component *)
f[x, y] > 0 && g[x, y] > 0 && y < 1 (* bottom component *)

Visualize:

cond = (f[x, y] < 0 && g[x, y] < 0 && y > 0) || (f[x, y] > 0 && g[x, y] > 0 && y < 1);

RegionPlot[cond, {x, -3, 3}, {y, -3, 3}]

enter image description here

We can ask for the area using an ImplicitRegion representation. Unfortunately (with Area) the exact form cannot be found, but we can get an approximate answer:

reg = ImplicitRegion[cond, {x, y}];

Area[reg, Method -> "NIntegrate"]
3.01605
Area[reg, Method -> {"NIntegrate", WorkingPrecision -> 20}]
3.0160539632460586365
$\endgroup$
  • $\begingroup$ That's just what I was looking for! I should have said that I was only taking the top region, so I just used NSolve to find the y value to split the condition, and got the top integral to 0.79... $\endgroup$ – Cameron Wood Nov 22 '17 at 18:06
  • $\begingroup$ Right. To use those endpoints you found with NSolve, then you'd need to do Integrate[1, {y, pt1, pt2}, {x, f1[y], f2[y]}], where f1 and f2 are messy inverses of f and g, respectively. $\endgroup$ – Chip Hurst Nov 22 '17 at 18:11
10
$\begingroup$

The area between two curves can be computed symbolically in a manner which makes it clear what is the basic idea behind the problem.

From the ContourPlot one can see that we could calculate x as a function of y for the both curves, we have:

Solve[y^3 - x^3 + 3 x^2 - 6 x - 4 y + 4 == 0, {y, x}, Reals] // Quiet
Solve[x^3 + y^3 + x - 2 y - 1 == 0, {y, x}, Reals] // Quiet
{{x -> Root[-4 + 4 y - y^3 + 6 #1 - 3 #1^2 + #1^3 &, 1]}}
  {{x -> Root[-1 - 2 y + y^3 + #1 + #1^3 &, 1]}}

Next we find the intersection points:

{r1, r2, r3} = y /. Solve[ Root[-4 + 4 y - y^3 + 6 #1 - 3 #1^2 + #1^3 &, 1]
                           == Root[-1 - 2 y + y^3 + #1 + #1^3 &, 1], y]
{ Root[67 - 382 #1 - 324 #1^2 - 95 #1^3 + 270 #1^4 + 216 #1^5 - 
          51 #1^6 - 72 #1^7 + 8 #1^9 &, 1], 
   Root[67 - 382 #1 - 324 #1^2 - 95 #1^3 + 270 #1^4 + 216 #1^5 - 
        51 #1^6 - 72 #1^7 + 8 #1^9 &, 2], 
   Root[67 - 382 #1 - 324 #1^2 - 95 #1^3 + 270 #1^4 + 216 #1^5 - 
        51 #1^6 - 72 #1^7 + 8 #1^9 &, 3]}

And here we have a symbolic solution:

Integrate[
  Abs[Root[-4 + 4 y - y^3 + 6 #1 - 3 #1^2 + #1^3 &, 1] - 
      Root[-1 - 2 y + y^3 + #1 + #1^3 &, 1]],
  {y, Root[67 - 382 #1 - 324 #1^2 - 95 #1^3 + 270 #1^4 + 216 #1^5 - 
           51 #1^6 - 72 #1^7 + 8 #1^9 &, 1], 
      Root[67 - 382 #1 - 324 #1^2 - 95 #1^3 + 270 #1^4 + 216 #1^5 - 
           51 #1^6 - 72 #1^7 + 8 #1^9 &, 3]}]

Since there we have roots of nineth order polynomial and the third order polynomila roots are quite involved the system cannot integrate the expression exactly. Nontheless we can find a good estimate:

NIntegrate[ Abs[Root[-4 + 4 y - y^3 + 6 #1 - 3 #1^2 + #1^3 &, 1] - 
                Root[-1 - 2 y + y^3 + #1 + #1^3 &, 1]], {y, r1, r3}]
 3.01605
Rotate [Plot[{ConditionalExpression[-Root[-4 + 4 y - y^3 + 6 #1 - 
    3 #1^2 + #1^3 &, 1], r1 <= y <= r3], 
             ConditionalExpression[-Root[-1 - 2 y + y^3 + #1 + #1^3 &, 1], 
                                    r1 <= y <= r3]},
          {y, -1.85, 1.55}, AxesLabel -> Automatic, 
          PlotRange -> {-1.85, 1.55}, AspectRatio -> Automatic, 
          Filling -> {1 -> {2}}], 90 Degree]

enter image description here

$\endgroup$
3
$\begingroup$

I propose another solution based on extracting the point of the mesh generated by the two contours. First, I generate the plots with the mesh curve $f=g=0$:

pl1 = Plot3D[#, {x, -2, 2}, {y, -2, 2}, Mesh -> {{0.}}, 
 MeshFunctions -> {#3 &}, PlotPoints -> 200, 
 MeshStyle -> {White, Thickness[0.007]}, ImageSize -> 500, 
 BoundaryStyle -> None, 
 PlotStyle -> {Red, Directive[Opacity[0.4]]}, 
 SphericalRegion -> True, AxesLabel -> Automatic] & /@ {f, g};

Now I extract the point of the both meshes and create a curve by interpolation:

ptsf = Join @@ Cases[Normal@#, Line[x1__] :> x1, Infinity] & /@ pl1;
flip1 = Interpolation[{#[[2]], #[[1]]} & /@ ptsf[[1]], 
Method -> "Spline", InterpolationOrder -> 3];
flip2 = Interpolation[{#[[2]], #[[1]]} & /@ ptsf[[2]], 
Method -> "Spline", InterpolationOrder -> 3];

The the region between the two curves (taking into account the intersection points):

Plot[{flip1[x], flip2[x]}, {x, solution[[1, 2]], solution[[3, 2]]}, 
Filling -> {1 -> {2}}, Frame -> True, FrameTicks -> True, 
RotateLabel -> False, 
FrameLabel -> (Style[#, 24, Italic,FontFamily -> "Times New Roman"] & /@ {"y", "x"}), 
PlotRange -> {{-2, 2}, All}]

enter image description here

It is easy obtain the area of the region for comparison to previous answers:

NIntegrate[Abs[flip1[x] - flip2[x]], {x, solution[[1, 2]], solution[[2, 2]], 
solution[[3, 2]]}, Method -> "GaussKronrodRule"]

yielding 3.01553. We can see that there is some difference in this value. I have checked that depends on the plot sampling, so coarse sampling gives lower value than those previous reported.

$\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.