Graphics and Plotting

Mathematica has a large number of built-in graphics programs for plotting and exploring mathematics graphically and for plotting and visualizing data.  Graphs and plots are effective ways of conveying complex information.  It is important to learn how to create and manipulate plots and graphics.  Mathematica is very rich in its graphical functions—we will explore just a small subset of all its capabilities.  We can get an idea of how many plotting routines are available by using a wildcard:

?*Plot*

While Mathematica has a large number of plotting routines, no one program does everything we need.  It can be useful to export numbers and graphics from Mathematica and operate on them individually or with other specialized programs.  The number of different creative graphical solutions available grows geometrically with the number of different graphical tools that are mastered.
Here, we will just produce a number of examples of plots and graphics...

Changing the Appearance of Simple Plots

Simple Plot with no options

Plot[Sin[x]/x,{x,-5 Pi,5Pi}]

Options[Plot]

PlotRange and PlotStyle

Plot[Sin[x]/x,{x,-5Pi,5Pi},PlotRange->{-0.25,1.25},PlotStyle->{Red,Thick}]

AxesLabel

Plot[Sin[x]/x,{x,-5Pi,5Pi},PlotRange->{-0.25,1.25},PlotStyle->{Red,Thick},AxesLabel->{"x","((Sin (x))/x)"}]

BaseStyle

Plot[Sin[x]/x,{x,-5Pi,5Pi},PlotRange->{-0.25,1.25},PlotStyle->{Red,Thick},AxesLabel->{"x","((Sin (x))/x)"},BaseStyle->{Large,FontFamily->"Helvetica",Italic}]

TicksStyle

Plot[Sin[x]/x,{x,-5Pi,5Pi},PlotRange->{-0.25,1.25},PlotStyle->{Red,Thick},AxesLabel->{"x","((Sin (x))/x)"},BaseStyle->{Large,FontFamily->"Helvetica",Italic},TicksStyle->{{Medium,Blue},{Medium,RGBColor[0.5,0.2,0]}}]

Controlling Precision and Size: The Mesh
(And an example of Interactive Control over Plot Appearance--Manipulate)

Mesh and MeshStyle

Plot[Sin[x]/x,{x,-5Pi,5Pi},PlotRange->All,PlotStyle->{Red,Thick},AxesLabel->{"x","((Sin (x))/x)"},BaseStyle->{Large,FontFamily->"Helvetica",Italic},TicksStyle->{{Medium,Blue},{Medium,RGBColor[0.5,0.2,0]}},Mesh->All,MeshStyle->{Black,PointSize[0.015]}]

MaxRecursion and PlotPoints

Plot[Sin[x]/x,{x,-5Pi,5Pi},PlotRange->All,PlotStyle->{Red,Thick},AxesLabel->{"x","((Sin (x))/x)"},BaseStyle->{Large,FontFamily->"Helvetica",Italic},TicksStyle->{{Medium,Blue},{Medium,RGBColor[0.5,0.2,0]}},Mesh->All,MeshStyle->{Black,PointSize[0.015]},MaxRecursion->2,PlotPoints->8]

Interactive Graphics: An Example of Manipulate

Manipulate[Plot[Sin[x]/x,{x,-5Pi,5Pi},PlotRange->All,PlotStyle->{Red,Thick},AxesLabel->{"x","((Sin (x))/x)"},BaseStyle->{Large,FontFamily->"Helvetica",Italic},TicksStyle->{{Medium,Blue},{Medium,RGBColor[0.5,0.2,0]}},Mesh->All,MeshStyle->{Black,PointSize[0.015]},MaxRecursion->recursion,PlotPoints->plotpoints],{{recursion,3},1,15,1},{{plotpoints,4},2,12,1}]

Plotting Multiple Curves and Plot Filling

Example of filling below a curve

Plot[Sin[x]/x,{x,-5Pi,5Pi},PlotRange->{-0.25,1.25},PlotStyle->{Red,Thick},AxesLabel->{"x","((Sin (x))/x)"},BaseStyle->{Large,FontFamily->"Helvetica",Italic},TicksStyle->{{Medium,Blue},{Medium,RGBColor[0.5,0.2,0]}},Filling->Automatic]

Combining several curves

Plot[{Sin[x]/x,Tan[x]/x},{x,-5Pi,5Pi},BaseStyle->{Thick}]

PlotStyles for each curve

Plot[{Sin[x]/x,Tan[x]/x},{x,-5Pi,5Pi},PlotStyle->{{Red,Thick},{Hue[0.3,1,.5],Thickness[0.005]}}]

Removing points with Exclusions

Plot[Tan[x]/x,{x,-5Pi,5Pi},BaseStyle->{Thick,Medium},Exclusions->{-Pi/2,Pi/2}]

Removing points with Table of Exclusions

Plot[Tan[x]/x,{x,-5Pi,5Pi},BaseStyle->{Thick,Medium},Exclusions->Table[p,{p,-9Pi/2,9Pi/2,Pi}]]

Multiple curves with exclusions

Plot[{Sin[x]/x,Tan[x]/x},{x,-5Pi,5Pi},PlotStyle->{{Red,Thick},{Hue[0.3,1,.5],Thickness[0.005]}},Exclusions->Table[p,{p,-9Pi/2,9Pi/2,Pi}]]

Filling  between curves

Plot[{Sin[x]/x,Tan[x]/x},{x,-5Pi,5Pi},PlotStyle->{{Red,Thick},{Hue[0.3,1,.5],Thickness[0.005]}},PlotRange->{-0.25,1.25},Exclusions->Table[p,{p,-9Pi/2,9Pi/2,Pi/2}],Filling->{2->{{1},{RGBColor[1,0,0,0.2],RGBColor[0,0,1,0.2]}}}]

Plotting Parametric Curves

ParametricPlot can produce a 2-D graph as a function of a single parameter, e.g. t. The cartesian coordinates, e.g. x and y, are specified as functions x(t) and y(t). Thus, a continuous variation of the single parameter t will trace out a trajectory in the x-y plane.  If two parameters are given, then ParametricPlot will fill the region x(r,t), y(r,t).

?ParametricPlot

Let's create a function that returns a list {x(t), y(t)} as a function of its argument t. We add a parameter n to make things interesting to look at...

MagicCircles[ t_, n_] :=  { Cos[n t -Pi +  2 Pi Quotient[n t,2Pi]/n ] + Cos[2 Pi Quotient[n t,2Pi]/n],  Sin[n t -Pi  +  2 Pi Quotient[n t,2Pi]/n ] + Sin[2 Pi Quotient[n t,2 Pi]/n]}

We use ParametricPlot on our function for n=5

ParametricPlot[MagicCircles[t,5],{t,0,2 Pi},PlotStyle->Thick,PlotRange->All]

We can use Manipulate to view how the plot changes with the parameter n, and also change the last point drawn to see how the curve gets traced out.

Manipulate[ParametricPlot[MagicCircles[t,ncirc],{t,0,lastp},PlotStyle->Thick,PlotPoints->6 ncirc,PlotRange-> 2.5{{-1,1},{-1,1}},Axes->False],{{ncirc,3},1,36,1},{{lastp,2 Pi},0.0001,2Pi}]

To demonstrate Parametrtic plot of a region, let's create a function that returns a list {x(r,t), y(r,t)} as a function of its arguments r and t. We add a parameter n to make things interesting to look at...

OrbitOrbit[ r_, t_, n_] :=  { r Cos[n t ] + Cos[t], r  Sin[n t] + Sin[t]}

If we fix r, this is a one dimensional curve as above.

ParametricPlot[Evaluate[OrbitOrbit[.5,t,12]],{t,-Pi, Pi},PlotStyle->Thick]

Manipulate[ParametricPlot[Evaluate[OrbitOrbit[.5,t,12]],{t,-Pi, endPoint},PlotRange-> 2{{-1,1},{-1,1}},PlotStyle->Thick],
{{endPoint,-3},-Pi,Pi}
]

Now we let both r and t vary. Some regions in the disk r ∈ (0.25,0.75) don't get covered, and others get covered one or more times.

ParametricPlot[Evaluate[OrbitOrbit[r,t,12]],{t,-Pi, Pi},{r,.25,.75},PlotStyle->{Thick,Red},Mesh->False,PlotPoints->72]

Just to demonstrate a pretty picture and to show  what kinds of things are possible without going into details, we use a Pure Function as the ColorFunction option.

ParametricPlot[Evaluate[OrbitOrbit[r,t,6]],{t,-Pi, Pi},{r,.25,.9},PlotStyle->{Thick,Red},Mesh->False,PlotPoints->36,ColorFunction->(Hue[#3,1,1,0.25]&)]

Plotting Data

ListPlot, PieChart, Histogram, Barchart, etc

The next command uses Mathematicas Integrated Data Resources, it will not retrieve the data unless you have an active internet connection

ElementData[]

Here is a list of properties that we can access from ElementData

ElementData["Properties"]

However, one should always question the provenence and accuracy of data... Let's make a sanity check: the stable phase of carbon at STP is graphite which is hexagonal (but not close packed).

ElementData[6,"StandardName"]
ElementData[6,"CrystalStructure"]

We create a list of the densities of the first one hundred elements. Data that is missing is reported with Missing[NotAvailable] or Missing[Unknown].

Densities = Table[ElementData[i,"Density"],{i,1,100}]

The data can be plotted with ListPlot

ListPlot[Densities]

ListPlot[Densities,BaseStyle->{Large,FontFamily->"Helvetica",PointSize[0.025]}]

To join the data, one can use either ListLinePlot

ListLinePlot[Densities,BaseStyle->{Large,FontFamily->"Helvetica",PointSize[0.025]}]
ListPlot[Densities,BaseStyle->{Large,FontFamily->"Helvetica",PointSize[0.025]}]

To see the data, we use the PlotMarkers Option.

ListLinePlot[Densities,BaseStyle->{Large,FontFamily->"Helvetica",PointSize[0.025]},PlotMarkers->Automatic,AxesLabel->{"Element Number","Density (MKS)"},ImageSize->Large]

Plotting Data with Screen Interaction

Example with Tooltip to make graphics interactive----put your mouse over a point and you get a pop-up with more information

ListLinePlot[Tooltip[Densities],BaseStyle->{Large,FontFamily->"Helvetica",PointSize[0.025]},PlotMarkers->Automatic,AxesLabel->{"Element Number","Density (MKS)"},ImageSize->Large]

This is a slightly more complicated example of Tooltip. We create a data structure with {x(i),y(i)} = {density(i), bulkmodulus(i)} and then tell Tooltip to pop-up the element's symbol when the mouse is over it.

ListPlot[Table[Tooltip[{ElementData[i,"Density"],ElementData[i,"BulkModulus"]},ElementData[i,"Abbreviation"],LabelStyle->{Large}],{i,1,100}],BaseStyle->{Large,FontFamily->"Helvetica",PointSize[0.025]},PlotMarkers->Automatic,AxesLabel->{"Density","Bulk Modulus"},PlotLabel->"MKS Units",ImageSize->Full]

Worked Example: Graphical  Data Exploring of the Elements

Here we do a small exercise to get a graphical representation of which Crystal Structures the elements form, and represent the frequency of each type. First we create a list of known elemental crystal structures for the first 100 elements.

CrystalStructures=Table[ElementData[i,"CrystalStructure"],{i,100}]

We manipulate our data to remove missing data (using Cases), and then use Tally to find how many times each Crystal structure appears.  We assign the results to a varialbe that we can use later to plot.

UniqueStructures = Tally[Cases[CrystalStructures,Except[Missing[_]]]]
MatrixForm[UniqueStructures]

Here is a bar chart showing the frequency of crystal structures. We had to manipulate the data (with Transpose) to get the numbers of each type and to get a list to hand to the Option BarLabels.

BarChart[Transpose[UniqueStructures][[2]],BarLabels->Transpose[UniqueStructures][[1]],BaseStyle->{Large,FontFamily->"Helvetica"},BarOrientation->Horizontal,ImageSize->Full]

We do the same thing, but with a pie chart.

PieChart[Transpose[UniqueStructures][[2]],PieLabels->Transpose[UniqueStructures][[1]],BaseStyle->{Large,FontFamily->"Helvetica"},ImageSize->Full]

As a last example, we produce a 3D histogram. The height of each bar corresponds to the number of elements in a range of melting points and range of densities.  We have to manipulate the data a bit to remove elements that have missing melting point OR density data.

histdata = DeleteCases[Table[{ElementData[i,"AbsoluteMeltingPoint"],ElementData[i,"Density"]},{i,100}],{Missing[_],_}|{_,Missing[_]}]
Histogram3D[histdata,{16,24},AxesLabel->{"Melting Point","Density","Number"} ,ImageSize-> Full]

Representations of 3D: Plot3D

This is a function to give the magnitude of the electrostatic potential at a point x,y,z for a point charge placed at xo,yo.

EPot[x_, y_ , z_ , xo_ , yo_] := (1/(Sqrt[(x-xo)^2 + (y-yo)^2 + z^2]))

SheetOLatticeCharge[x_, y_ , z_] := Sum[EPot[x,y,z,xo,yo],{xo,-5,5},{yo,-5,5}]

SheetOLatticeCharge represents the electric field produced by an 11 by 11 array of point charges arranged on the x-y plane at z = 0. The following command evaluates and plots the field variation in the plane z = 0.25:

Plot3D[Evaluate[SheetOLatticeCharge[x,y,0.25]],{x,-6,6},{y,-6,6}]

Note below how theplot is set to contain the output of the Plot3D command---it is now a symbol assigned to a graphics object.  The number of plotpoints is increased so that we can resolve all the bumps.  This will take a while to compute on most machines.

theplot =Plot3D[Evaluate[SheetOLatticeCharge[x,y,0.25]],{x,-6,6},{y,-6,6},PlotPoints->60]

This demonstrates the use of RegionFunction plot option which is pure function. Here, only the region inside a cylinder with radius 9 ((x^2) + (y^2) ≤ (9^2)) is plotted.

Plot3D[Evaluate[SheetOLatticeCharge[x,y,0.25]],{x,-9,9},{y,-9,9},PlotPoints->60,RegionFunction->(#1^2 + #2^2 ≤ 81&)]

This demonstrates the use of the ColorFunction  plot option which is pure function. Here we use one of Mathematica ColorData functions.

Plot3D[Evaluate[SheetOLatticeCharge[x,y,0.25]],{x,-9,9},{y,-9,9},PlotPoints->60,RegionFunction->(#1^2 + #2^2 ≤ 81&),ColorFunction->(ColorData["TemperatureMap"][#3]&)]

Lower Dimensional Representations of 3D: ContourPlot

theconplot =ContourPlot[Evaluate[SheetOLatticeCharge[x,y,0.25]],{x,-6,6},{y,-6,6},PlotPoints->32]

theconplot =ContourPlot[Evaluate[SheetOLatticeCharge[x,y,0.25]],{x,-4,4},{y,-4,4},PlotPoints->50,ColorFunction->Hue,Contours->24]

thedenplot =DensityPlot[Evaluate[SheetOLatticeCharge[x,y,0.25]],{x,-4,4},{y,-4,4},PlotPoints->50,ColorFunction->ColorData["GreenBrownTerrain"]]

Graphics Primitives

It can be useful to be able to build up arbitrary graphics objects piece-by-piece using simple "graphics primitives" like Circle:

thecirc =Graphics[Circle[{2,2},1.5]]

Show[thecirc,Axes->True]

Show[thecirc,Axes->True,AxesOrigin->{0,0},AspectRatio->1]

Now we take a simple plot…

cosplot = Plot[Cos[x],{x,0,4Pi}]

Adding Graphics Primitives to Plots (or other graphics objects) using the built-in Drawing Tool

Mathematica6 now has a simple drawing editor that allows you add text, arrows, lines, and shapes to existing graphics. To do this, select the previous graphics output for the cosine plot.  While the graphics are selected, use the Menu Item "Drawing Tools" under Graphics. After you have added shapes, text, etc.. move the cursor to the left of the selected graphics object and type a symbol (below, I used "thenewplot") for the new (combined) graphics object to be assigned to.

thenewplot=

thenewplot

Combining Graphical Objects using Show.

and overlay some text in places of our own choosing…

Show[cosplot,Graphics[Text["One Wavelength",{2 Pi,0.5}]],Graphics[Text["Two\nWavelengths",{4 Pi,0.5}]],PlotRange->All]

Show[thenewplot,Graphics[Text["One Wavelength",{2 Pi,1.1}]],Graphics[Text["Two Wavelengths",{4 Pi,1.1}]],PlotRange->All]

A Worked Example using Graphics Primitives: The Wulff Construction

This next example shows a clever way to perform a famous thermodynamic graphical construction called  the  Wulff construction.  For surface tensions that are a function of their orientation, the Wulff construction gives the enclosing surface of lowest total energy.  The following code creates the graphics objects corresponding to the Wulff construction; it will be combined with an orientation-dependent surface tension to produce a prediction of the energy-minimizing shape

wulffline[{x_,y_},wulfflength_] := Module[{θ,wulffhalf=wulfflength*0.5,x1,x2,y1,y2},θ=ArcTan[x,y];
 x1 = x + wulffhalf*Cos[θ + Pi/2];x2=  x + wulffhalf*Cos[θ - Pi/2];y1 = y + wulffhalf*Sin[θ + Pi/2];y2= y + wulffhalf*Sin[θ - Pi/2];
Graphics[Line[{{x1,y1},{x2,y2}}]]
]

These last two plots have to do with the orientation dependence of the surface energy of a 2-D crystal (GammaPlot)and the equilibrium shape that a single crystal having this surface energy will adopt to minimize its total surface energy.

gammaplot[ theta_ , anisotropy_ , nfold_] := {Cos[theta]+ anisotropy*Cos[(nfold+1)*theta],Sin[theta] + anisotropy*Sin[(nfold+1)*theta]}

GammaPlot = ParametricPlot[gammaplot[t,0.1,4],{t,0,2 Pi},PlotStyle->{{Thickness[0.01],RGBColor[1,0,0]}}]

Show[Table[wulffline[gammaplot[t,0.1,4],2],{t,0, 2 Pi, 2 Pi/100}],GammaPlot]

The inner envelope of the tangent lines to the GammaPlot in the figure above illustrates the energy-minimizing shape of the crystal. For this particular GammaPlot, the shape is approximately a square but the corners are rounded.

ToutesDesLoups[anisotropy_, nfold_] := Module[{GammaPlot },GammaPlot= ParametricPlot[gammaplot[t,anisotropy,nfold],{t,0,2 Pi},PlotStyle->{{Thickness[0.01],RGBColor[1,0,0]}}];
Show[Table[wulffline[gammaplot[t,anisotropy,nfold],3],{t,0, 2 Pi, 2 Pi/100}],GammaPlot]]

Manipulate[ToutesDesLoups[aniso,nfold],{{aniso,0.1},-0.9,0.9},{{nfold,6},2,16,1}]

Animations

Producing animations are very simple. We have already met some examples using Manipulate which provides lower-level control over animations, but anything that produces can be turned into an animation by animating over a parameter. Here are some simple examples:

Fxt[x_, t_] := Sin[3 (x+10 -t)]Exp[-(x + 10 - t)^2] - Sin[3(x -10 + t)] Exp[-(x - 10 + t)^2]

Here is a simple plot at "time" = 0

Plot[Fxt[xvar,0],{xvar,-15,15},PlotRange->{-1,1},PlotStyle->{Thick,Red},Filling->Axis,FillingStyle->{RGBColor[0,0.5,0,0.5],RGBColor[0,0,0.5,0.5]}]

To Animate :

Animate[Plot[Fxt[xvar,timevar],{xvar,-15,15},PlotRange->{-1,1},PlotStyle->{Thick,Red},Filling->Axis,FillingStyle->{RGBColor[0,0.5,0,0.5],RGBColor[0,0,0.5,0.5]}],{timevar,0,25}]

This is the solution to the temperature evolution equation (the diffusion equation) for a square of length L initially at 500K embedded in a plate initially at 100K , κ is the themal diffusivity (units (length^2)/time).  We introduce a "normalized" time and space variables variable τ= κ t/(L^2) and ξ = x/L and η=y/L

TempSquare = 100 +  400 Integrate[((Exp[-(((((x-xo))^2) + (((y-yo))^2))/(4 κ t))])/(4π κ t)),{xo,-L/2,L/2},{yo,-L/2,L/2}]
NormalizeRules ={t -> τ (L^2)/κ,  x ->ξ L ,  y-> η L ,  xo-> ξo L,  yo-> ηo L};
TempSquare =Simplify[TempSquare/.NormalizeRules,Assumptions->κ>0 && L>0]

We divide by 500 so that the temperatures should scale between zero and one, and then use ColorFunctionScaling->False so that the colors are consistent over time.

ListAnimate[Table[Plot3D[TempSquare/500,{η,-1,1}, {ξ,-1,1},PlotRange->{0,1},PlotPoints->50,ColorFunction->"TemperatureMap",ColorFunctionScaling->False],{τ,0.001, .1,0.002}]]

Worked Animation Example with Graphics Primitives: Visualizing the Random Walk

Many processes in the natural world can be modelled by random walks. An example is the trajectory of a particular  molecule in a gas at low pressure.  The next few lines of code produce random walks and display them graphically:

randomwalk[0] = {0,{0,0}}

randomwalk[0][[2]]

randomwalk[nstep_Integer?Positive]:=randomwalk[nstep]={nstep,randomwalk[nstep-1][[2]]+RandomReal[0.5] {Cos[theta= RandomReal[2 π]],Sin[theta]}}

randomwalk[4]

Create a function that returns a graphic object putting the step number at the correct place:

gtext[nstep_Integer?NonNegative] := gtext[nstep] = Graphics[Text[ToString[randomwalk[nstep][[1]]],randomwalk[nstep][[2]]]];

locations=Show[Table[gtext[i],{i,0,100}],PlotRange->All,AspectRatio->1]

Create a function that returns a graphic object connecting a line between the current position and the previous.

gline[nstep_Integer] := gline[nstep] = Graphics[Line[{randomwalk[nstep-1][[2]],randomwalk[nstep][[2]]}]];

jumps =Show[Table[gtext[i],{i,0,100}],Table[gline[j],{j,1,100}],PlotRange->All,AspectRatio->1]

We can animate the jumps, but the result will not be satisfactory because each frame will be a different size...

Animate[Show[gtext[i],gline[i]],{i,1,49,1}]

If we use the PlotRange from a graphical object that contains all the points, we can fix the framesize, we use AbsoluteOptions

prange= PlotRange/.AbsoluteOptions[locations]

Animate[Show[gtext[i],gline[i],PlotRange->prange],{i,1,100,1}]

The animation is more satisfactory if we include the history... here we create a table of objects leading up to the current step.

Animate[Show[Table[{gtext[i],gline[i]},{i,1,j}],PlotRange->prange],{j,2,100}]

A Worked Example: Common Tangent Construction and Spinodal Points for a Single Curve 
(Part A: The Spinodal Region)

A prototype molar free energy of mixing using the same xlogx function for the ideal entropy of mixing terms.  The temperature term is a scaled energy (RT), and it is assumed that enthalpies have been scaled so that the temperatures of interest (if there are any) are between T=0 and T=10.

xlogx[0]=xlogx[1]=xlogx[0.0] = xlogx[1.0]=0;xlogx[x_] := x Log[x]
Gmolar[X_, T_] := 5 X (1-X) + T (xlogx[X] +  xlogx[1-X]) + X/2

Here is the shape of our prototype free energy at T=3/2

p1 =Plot[Gmolar[x,3/2],{x,0,1},PlotStyle->Thick]

p1 =Plot[Gmolar[x,3/2],{x,0,1},PlotStyle->Thick,BaseStyle->{Large},AxesLabel-> {"Molar Composition","(OverBar[G])"}]

We will need the bounds of the above graphics object:

{{graphxmin,graphxmax},{graphymin,graphymax}}= PlotRange/.AbsoluteOptions[p1,PlotRange]

First let's determine where the spinodal region (by finding where the second derivative with respect to composition is negative

ddg =D[Gmolar[x,3/2],{x,2}]

Then, use RegionPlot to illustrate the range over which spinodal decomposition is spontaneous

p2 =RegionPlot[ddg<0, {x,graphxmin,graphxmax},{T,graphymin,graphymax},PlotStyle->RGBColor[0,1,.5,0.1]]

Show them both together to identify the spinodal region

Show[p1,p2]

A Worked Example: Common Tangent Construction and Spinodal Points  (Part B: Computing the Common Tangent)

We can use the ConvexHull to find the common tangent lines; this function is in the Computational Geometry Package.

<<ComputationalGeometry`

First we compute a list of values along the molar free energy curve, then compute those that lie outside the common tangent(s) (i.e., the convex hull). Because the points are given in order, we might as well
sort them on the way back out.  Note, the convex hull program gives the indices of the vertices that are on the hull.

npoints = 100;
gvals =  Table[{x,Gmolar[x,3/2]},{x,0,1,1/N[npoints-1]}];

We only want the lower convex hull; therefore we add some "fictive" points to the beginning and the end of the data.  The the fictive points add a rectangle to the top of the curve that should be part of the computed convex hull.

gmax = Max[Transpose[gvals][[2]]];
PrependTo[gvals,{0,10*Abs[gmax]}];
AppendTo[gvals,{1,10*Abs[gmax]}];

After we compute this hull, we shift the hull by one and take off its first and last element. We strip the first and last element from the discrete values of free energy as well.

chull = Sort[ConvexHull[gvals]];
chull = Drop[Drop[chull -1,1],-1]
gvals = Drop[Drop[gvals,1],-1]

The common tangent(s) correspond to gaps in the vertex list of the common tangent. We will use Split to find the set of continous sequences.

convexparts = Split[chull,(#2 - #1<2)&]

A Worked Example: Common Tangent Construction and Spinodal Points (Part C: Creating Graphics Objects and Combining Plots)

Now we create graphics objects for each of the two-phase regions (i.e., the gaps in the convex hull) and collect them all into a graphics list for subsequent display.

len=Length[convexparts];
graphicslist={};
i=1;
While[i+1 ≤len,leftpoint = gvals[[Last[ convexparts[[i]] ] ]] ;
 rightpoint =gvals[[ First[ convexparts[[i+1]] ] ]];
ctline = {Red, Thick, Line[{leftpoint, rightpoint}]};
twophaseregion = {RGBColor[0.5,0,0,0.2],Rectangle[{leftpoint[[1]],graphymin},{rightpoint[[1]],graphymax}]};
AppendTo[graphicslist,ctline];
AppendTo[graphicslist,twophaseregion];
i++
]
p3=Graphics[graphicslist]

Show[p1,p2,p3]

A Worked Example: Common Tangent Construction and Spinodal Points  (Part D: Wrapping it up in a Tidy Module)

Put all the steps together in a Module for later use. Let the module take a function of two arguments for the molar free energy as a function of X and T, and an optional third parameter which specifies the number of data points to compute on the hull.

Needs["ComputationalGeometry`"];

CommonTangentConstruction[Gm_,T_, npts_:100]:=
Module[{x,y,p1,p2,p3, gxmin,gxmax,gymin,gymax,ddg,gvals,gmax,chull,conprts,len,glist={},i=1,lftpt,rtpt,ctline,twophasreg}, 
p1 =Plot[Gm[x,T],{x,0,1},PlotStyle->Thick];
{{gxmin,gxmax},{gymin,gymax}}= PlotRange/.AbsoluteOptions[p1,PlotRange];
ddg =D[Gm[x,T],{x,2}];
p2 =RegionPlot[ddg<0, {x,gxmin,gxmax},{y,gymin,gymax},PlotStyle->RGBColor[0,1,.5,0.1]];
gvals =  Table[{x,Gm[x,T]},{x,0,1,1/N[npts-1]}];
gmax = Max[Transpose[gvals][[2]]];
PrependTo[gvals,{0,10*Abs[gmax]}];
AppendTo[gvals,{1,10*Abs[gmax]}];
chull = Sort[ConvexHull[gvals]];
chull = Drop[Drop[chull -1,1],-1];
gvals = Drop[Drop[gvals,1],-1];
conprts = Split[chull,(#2 - #1<2)&];
len=Length[conprts];
While[i+1 ≤len,lftpt = gvals[[Last[ conprts[[i]] ] ]] ;
 rtpt =gvals[[ First[ conprts[[i+1]] ] ]];
ctline = {Red, Thick, Line[{lftpt, rtpt}]};
twophasreg = {RGBColor[0.5,0,0,0.2],Rectangle[{lftpt[[1]],gymin},{rtpt[[1]],gymax}]};
AppendTo[glist,ctline];
AppendTo[glist,twophasreg];
i++
];
p3=Graphics[glist];
Show[p1,p2,p3]
]

A Worked Example: Common Tangent Construction and Spinodal Points  (Part E: Visualize the Fruit)

CommonTangentConstruction[Gmolar,1.5]

Manipulate[CommonTangentConstruction[Gmolar,T,300],{{T,2},0,3}]