Plotting over non-rectangular regions

Initialization

<< Graphics`

\(plot1\  = 
      FilledPlot[{4\  - \ x\^2, 1}, {x, \(-\@3\), \@3}, 
        PlotRange\  -> \ {{\(-\@3\), 2.3}, {0, 4}}, 
        AspectRatio\  -> \ Automatic, 
        Epilog\  -> \ {Text["y = 1", { .6,  .8}], 
            Text[\*""y = 4 - (x^2)"", {1.7, 3.5}]}, 
        PlotLabel\  -> \ "Region Determined by Two
  Intersecting Curves", Ticks\  -> \ {Automatic, {1, 2, 3, 4}}, 
        DisplayFunction\  -> \ Identity]\)

\(plot2\  = 
      FilledPlot[{4\  - \ x\^2, Abs[x]/\@3}, {x, \(-\@3\), \@3}, 
        PlotRange\  -> \ {{\(-2\), 2.3}, {0, 4}}, 
        AspectRatio\  -> \ Automatic, 
        Epilog\  -> \ {Text[\*""y = x/(Sqrt[3])"", {1.3,  .3}], 
            Text[\*""y = -x/(Sqrt[3])"", {\(-1.4\),  .3}], 
            Text[\*""y = 4 - (x^2)"", {1.7, 3.5}]}, 
        PlotLabel\  -> \ "Region Determined by a Curve
  and Two Straight Lines", Ticks\  -> \ {Automatic, {1, 2, 3, 4}}, 
        DisplayFunction\  -> \ Identity]\)

botfn[x_] := Which[x ≤  -(Sqrt[2]),x/((Sqrt[2]) - (Sqrt[3])) +  (Sqrt[2])/((Sqrt[2]) - (Sqrt[3])),-(Sqrt[2]) < x < (Sqrt[2]),(Sqrt[2 - (x^2)]),True,-x/((Sqrt[2]) - (Sqrt[3]))+ (Sqrt[2])/((Sqrt[2]) - (Sqrt[3]))]

\(plot3\  = 
      FilledPlot[{botfn[x], 4\  - \ x\^2}, {x, \(-\@3\), \@3}, 
        PlotRange\  -> \ {{\(-2\), 2.3}, {0, 4}}, 
        AspectRatio\  -> \ Automatic, 
        Epilog\  -> \ {Text[\*""y = \!\(\@\(2\\  - \\ x\^2\)\)"", { .82, 1.6}], 
            Text[\*""y = 4 - (x^2)"", {1.7, 3.5}]}, 
        PlotLabel\  -> \ "Region Determined by two Curves
    and Two Straight Lines", Ticks\  -> \ {Automatic, {1, 2, 3, 4}}, 
        DisplayFunction\  -> \ Identity]\)

xvar[a_, b_]  :=  -1.4142135623730951` - 
        0.31783724519578205` a + (2.8284271247461903` + 0.6356744903915641` a) b

yvar[a_, 
        b_] :=  -2.8284271247461903` (-1 + 
              a) (Sqrt[-(-1.` + b) (-5.5511151231257815`*^-17 + b)]) - 
        11.999999999999998` a (-1.0773502691896257` + 
              b) (0.07735026918962581` + 
              b)

crvtbl  = 
        Chop[ Table[{xvar[a, b], yvar[a, b]}, {a, 0, 1,  .1}, {b, 0, 
              1,  .04}], (10^(-7))];

lntbl  =  
        Chop[Table[{xvar[a, b], yvar[a, b]}, {b, 0, 1,  .04}, {a, 0, 
              1,  .1}], (10^(-7))];

An overview

Plotting over non-rectangular regions can always be accomplished using ParametricPlot3D, but it is often difficult to create the plot without some involved calculations. RegionPlot3D does the work for you if you can define the boundaries of the region by functions of the form y = f(x) or x = g(y).

How the function works

Suppose we have a region defined by the two curves y = 4 - (x^2)for -(Sqrt[3])<  x < (Sqrt[3]) and y = (Sqrt[2  -  (x^2)])for -(Sqrt[2])<  x < (Sqrt[2])and the straight lines joining the ends of the curves. (See Figure 1.)

\(FilledPlot[{botfn[x], 4\  - \ x\^2}, {x, \(-\@3\), \@3}, 
      PlotRange\  -> \ {{\(-2\), 2.3}, {0, 4}}, 
      AspectRatio\  -> \ Automatic, 
      Epilog\  -> \ {Text[\*""y = \!\(\@\(2\\  - \\ x\^2\)\)"", { .64,  .45}], Text[\*""y = 4 - (x^2)"", {1.4, 3.5}]}, 
      PlotLabel\  -> \ "Figure 1", 
      Ticks\  -> \ {Automatic, {1, 2, 3, 4}}, 
      ImageSize\  -> \ 280]\)

We proceed as follows. First we divide the x-range both curves into 20 equal parts. The divisions will be of different size for the two curves. Next we mark off the points on the curves determined by these subdivisions and join corresponding points on the two curves with straight lines. Next we divide the length of each of the straight lines into ten equal units. This creates a grid which is shown in Figure 2. Then all we need do to create a 3D plot of a surface over this region is to assign a third coordinate to each of the grid points. This is what RegionPlot3D does.

\(Show[
      Graphics[{{Hue[ .95], Line /@ Take[crvtbl, {2, 11}], 
            AbsoluteThickness[2], 
            Line /@ {First[crvtbl], Last[crvtbl]}}, {Hue[ .7], 
            Line /@ lntbl}}], AspectRatio\  -> \ Automatic, 
      Axes\  -> \ True, PlotLabel\  -> \ "Figure 2", 
      Epilog\  -> \ {Text[\*""y = \!\(\@\(2\\  - \\ x\^2\)\)"", { .6,  .7}], Text[\*""y = 4 - (x^2)"", {1.25, 3.5}]}, 
      ImageSize\  -> 260]\)

An example

Here is an example plotting the function z = y/3.   - (x^2)/5.  over the above region.

RegionPlot3D[{fn_, {var1_, var2_}}, {f1_, {var_, varmin1_, 
            varmax1_}}, {f2_, {var_, varmin2_, varmax2_}}, opts___] := 
      Module[{scalefactor1, scalefactor2, myvar1, myvar2, tt, ffn, ff1, ff2, 
          rr, pltopts, styleopts, t, r}, 
pltopts  =  
          Select[{opts}, #[[1]] =!= 
                SurfaceStyle &]; 
styleopts  =  SurfaceStyle /. {opts} /. 
            SurfaceStyle  ->  {}; 
ffn[tt_, rr_] = 
          If[Head[fn] === Symbol && ! MemberQ[{var1, var2, var}, fn], 
            fn[tt, rr], 
            fn /. {var1 -> tt, 
                var2 ->  rr}]; 
ff1[tt_] = 
          If[Head[f1] === Symbol && ! MemberQ[{var1, var2, var}, fn1], 
            f1[tt], f1 /. var -> tt]; 
ff2[tt_] = 
          If[Head[f2] === Symbol && ! MemberQ[{var1, var2, var}, fn2], 
            f2[tt], f2 /. var -> tt]; 
Which[
          var === var1, 
scalefactor1[
              r_]  =  (varmax1 - varmin1) r  + 
              varmin1; 
scalefactor2[
              r_]  =  (varmax2 - varmin2) r  +  
              varmin2; 
myvar1[t_, r_]  =  
            Simplify[(1 - t) scalefactor1[r]  + 
                t scalefactor2[r]]; 
myvar2[t_, r_]  =  
            Simplify[(1 - t) ff1[scalefactor1[r]]  +  
                t ff2[scalefactor2[r]]]; 
ff2n[tt_, rr_] = 
            Simplify[
              ffn[myvar1[tt, rr], myvar2[tt, rr]]]; 
pltfn[
              t_, r_] = {myvar1[t, r], myvar2[t, r], ff2n[t, r], 
              styleopts}, 
var  ===  
            var2, 
scalefactor1[
              r_]  =  (varmax1 - varmin1) r  + 
              varmin1; 
scalefactor2[
              r_]  =  (varmax2 - varmin2) r  +  
              varmin2; 
myvar2[t_, r_]  =  
            Simplify[(1 - t) scalefactor1[r]  + 
                t scalefactor2[r]]; 
myvar1[t_, r_]  =  
            Simplify[(1  -  t) ff1[scalefactor1[r]]  +  
                t ff2[scalefactor2[r]]]; 
ff2n[tt_, rr_] = 
            Simplify[
              ffn[myvar1[tt, rr], myvar2[tt, rr]]]; 
pltfn[
              t_, r_] = {myvar1[t, r], myvar2[t, r], ff2n[t, r], 
              styleopts}, 
True, 
scalefactor1[t_]  =  (varmax1 - varmin1) t  + 
              varmin1; 
scalefactor2[
              t_]  =  (varmax2 - varmin2) t  +  
              varmin2; 
basecoords[r_, t_]  =  
            Simplify[(1 - r) ff1[scalefactor1[t]] +  
                r ff2[scalefactor2[t]]]; 
pltfn[r_, 
              t_]  =  
            Flatten[{basecoords[r, t], 
                Simplify[Apply[ffn, basecoords[r, t]]], 
                styleopts}]
]; 
ParametricPlot3D[
          Chop[Evaluate[pltfn[r, t]], 10^-7], {t, 0, 1}, {r, 0, 1}, 
          Evaluate[pltopts], 
          Compiled  ->  False]
]

RegionPlot3D[{y/3.   - (x^2)/5. , {x, 
          y}}, {(Sqrt[2.   -  (x^2)]), {x, -(Sqrt[2]), (Sqrt[2])}}, {4  -  
          (x^2), {x, -(Sqrt[3]), (Sqrt[3])}}, ImageSize -> 260, 
      PlotPoints  ->  {15, 45}]

Using the function

Boundaries of the form y = f(x)

The input for this function follows the basic form for all Mathematica plot functions. It is always possible to combine several plots and so we need only work with three types of regions. These are illustrated in the Figure 3. The function will not work correctly if any of the interior angles are greater than π radians (180 degrees).

\(Show[
      GraphicsArray[{{plot1, plot2}, {plot3}}, Frame\  -> \ True, 
        PlotLabel\  -> \ 
          StyleForm["         Figure 2
Possible plot regions", "Subsection"]], ImageSize\  -> \ 520, 
      DisplayFunction\  -> \ $DisplayFunction]\)

Arguments for RegionPlot3D

The arguments for RegionPlot3D for this type of region are (the variable specification in the second and third arguments must be the same.):

RegionPlot3D[{surface function, {var1, 
            var2}}, {boundary function 1, {var, varmin1, 
            varmax1}}, {boundary function2, {var, varmin2, varmax2}}, 
        options]

Case one

The three region examples given above all share the property that the boundaries can be written as functions of x. We start with the first case by finding the intersection of the two curves. You can often do this by inspection but you can also do it using Mathematica.

sol1  =  Solve[4  -  (x^2) ==  1, x]

pts  =  Map[x /. # &, sol1]

We use Flatten on the second and third arguments to get them in the proper form.

RegionPlot3D[{ y/2, {x, y}}, {4.  - (x^2), Flatten[{x, pts}]}, {1. , 
        Flatten[{x, pts}]}]

Case Two

Here we take one of the curves to be a point. Everything else is handled in the same way as in the first example.

RegionPlot3D[{y, {x, y}}, {0, {x, 0, 0}}, {4.  - (x^2), 
        Flatten[{x, pts}]}]

Case Three

Here things are handled as in the first case except the end points are different for each of the boundaries. Since the boundary curves are given as functions of x we give the points in ascending order of x-coordinates.

RegionPlot3D[{y/3.   - (x^2)/5. , {x, 
          y}}, {(Sqrt[2.   -  (x^2)]), {x, -(Sqrt[2]), (Sqrt[2])}}, {4  -  
          (x^2), {x, -(Sqrt[3]), (Sqrt[3])}}, ImageSize -> 260, 
      PlotPoints  ->  {15, 15}]

In this case we can also use the case one situation by using Which to define the lower boundary of the region. Notice how the surface mesh changes when we do this.

botfn[x_] := Which[x ≤  -(Sqrt[2.]) ,x/((Sqrt[2.])  - (Sqrt[3.]) ) +  (Sqrt[2.]) /((Sqrt[2.])  - (Sqrt[3.]) ),-(Sqrt[2.])  < x < (Sqrt[2.]) ,(Sqrt[2.  - (x^2)]),True,-x/((Sqrt[2.])  - (Sqrt[3.]) )+ (Sqrt[2.]) /((Sqrt[2.])  - (Sqrt[3.]) )]

RegionPlot3D[{y/3.   - (x^2)/5. , {x, y}}, {botfn[
          x], {x, -(Sqrt[3]), (Sqrt[3])}}, {4  -  (x^2), {x, -(Sqrt[3]), (Sqrt[3])}}, 
      ImageSize -> 260, PlotPoints  ->  15]

Options

RegionPlot3D will take all of the standard Options of ParametricPlot3D. It also takes a new Option SurfaceStyle. SurfaceStyle takes any value which is a valid argument for the fourth coordinate in ParametricPlot3D. Thus to eliminate the "grid" on our surface we could use SurfaceStyle ->  EdgeForm[]. We can also modify the color of the surface by adding color directives. SurfaceColor -> {SurfaceColor[RGBColor[.45,.45,0],RGBColor[.2,.2,.4],10],EdgeForm[]}.

plot1  = 
      RegionPlot3D[{y/3.   - (x^2)/5. , {x, 
            y}}, {(Sqrt[2.   -  (x^2)]), {x, -(Sqrt[2]), (Sqrt[2])}}, {4  -  
            (x^2), {x, -(Sqrt[3]), (Sqrt[3])}}, ImageSize -> 260, 
        PlotPoints  ->  {15, 45}, 
        SurfaceStyle  ->  {SurfaceColor[RGBColor[ .45,  .45, 0], 
              RGBColor[ .2,  .2,  .4], 10], EdgeForm[]}, 
        ViewPoint -> {-2.606,  -1.381,  1.659}, 
        LightSources -> {{{1.`, 0.`, 1.`}, 
              RGBColor[1.`, 1. , 1. ]}, {{1.`, 1.`, 1.`}, 
              RGBColor[1. , 1. , 1.`]}, {{0.`, -1.`, 1.`}, 
              RGBColor[1. , 1.`, 1.`]}}]

RegionPlot3D also works in conjunction with other Graphic functions such as shadow.

Shadow[plot1, XShadow  ->  False, YShadow  ->  False, 
      ViewPoint -> {-2.606,  -1.381,  1.659}, 
      ZShadowPosition  ->  -7/4, Boxed  ->  False, 
      Axes  ->  True, 
      LightSources -> {{{1.`, 0.`, 1.`}, 
            RGBColor[1.`, 1. , 1. ]}, {{1.`, 1.`, 1.`}, 
            RGBColor[1. , 1. , 1.`]}, {{0.`, -1.`, 1.`}, 
            RGBColor[1. , 1.`, 1.`]}}]

Boundaries of the form x  = g(y)

Alternate Arguments for RegionPlot3D

In some cases it is easier to express the boundaries of the region as a function of y. RegionPlot3D accepts this configuration and automatically orients the plot correctly. Here is an example creating a section of a punctured sphere.

Quarter sphere

ff1[y_] := (Sqrt[1 - (y^2)])

ff2[y_]  :=  (Sqrt[4  -  (y^2)])

ff4[x_, y_]  := (Sqrt[4 - ((x^2)  +  (y^2))])

plot2  = 
      RegionPlot3D[{ff4, {x, y}}, {ff1, {y, -1, 1}}, {ff2, {y, -2, 
            2}}, PlotPoints  ->  {30, 35}, 
        ViewPoint -> {-1.812,  2.548,  1.294}, 
        SurfaceStyle  ->  EdgeForm[], 
        AxesLabel  ->  {"x-axis", "y-axis", ""}]

We can get the other half of the above sphere by defining a second set of boundary functions.

ff1N[y_] := -(Sqrt[1 - (y^2)])

ff2N[y_]  :=  -(Sqrt[4  -  (y^2)])

plot3  = 
      RegionPlot3D[{ff4, {x, y}}, {ff1N, {y, -1, 1}}, {ff2N, {y, -2, 
            2}}, PlotPoints  ->  {30, 35}, 
        ViewPoint -> {-1.812,  2.548,  1.294}, 
        SurfaceStyle  ->  EdgeForm[], 
        AxesLabel  ->  {"x-axis", "y-axis", ""}]

Combining plots

And here we put them together.

Show[{plot2, plot3}]

By including a cylinder of appropriate size we can make the sphere appear to be solid. We use the following to determine the height of the cylinder.

ff4[1, 0]

plot4  =  
      ParametricPlot3D[{Cos[t], Sin[t], z, EdgeForm[]}, {t, 0, 2 π}, {z, 
          0, (Sqrt[3])}, PlotPoints  ->  {36, 4}, 
        DisplayFunction  ->  Identity]

Show[{plot2, plot3, plot4}, Boxed  ->  False]

Boundaries in polar coordinates

The arguments for RegionPlot3D for this type of region are :

RegionPlot3D[{surface function, {var1, var2}}, {{x1[var], 
            y1[var]}, {var, varmin1, varmax1}}, {{x2[var], y2[var]}, {var, 
            varmin2, varmax2}}, options]

In the above example of the sphere it would be easier to give the boudary curves in polar coordinates. RegionPlot3D will also accomodate this case, but the form of the arguments changes slightly. The boundary functions must be given as ordered pairs of functions and the variable name cannot be either one of the variables given in the first argument. The first function used to define the surface should still be given in cartesian coordinates.

plot5  =  
      RegionPlot3D[{ff4, {x, y}}, {{Cos[t], Sin[t]}, {t, 0, 
            2 π}}, {{2 Cos[t], 2 Sin[t]}, {t, 0, 2 π}}, 
        PlotPoints  ->  {30, 35}, 
        ViewPoint -> {-1.812,  2.548,  1.294}, 
        SurfaceStyle  ->  EdgeForm[], 
        AxesLabel  ->  {"x-axis", "y-axis", ""}]

Show[{plot4, plot5}, 
      DisplayFunction  ->  $DisplayFunction]

Limitations of the function

There are some limitations of the function due to the accumulation of numerical error. First the independent variable used to define the boundaries is scaled. The scaled variable is then fed to boundary functions. Both of these functions are then used as arguments in the surface function leading to the following composition of functions:

	{scalefunc[x], boundaryfunc[scalefunc[x]], surfacefunc[scalefunc[x],boundaryfunc[scalefunc[x]]]}
			
Since all of the plot functions in Mathematica automatically use machine precision numbers it is inevitable that with some boundary functions small imaginary parts are going to arise. I have attempted to account for this by setting the critical value for Chop to (10^(-7))which works in most cases. This may, however, lead to problems with small regions where more decimal places are needed to define the boundaries of the region causing necessary values to be discarded.