The PrecisionPlot package
by Ted Ersek  (ErsekTR@navair.navy.mil)

Notebook and package revised on  30 March 1999

Motivation for defining PrecisionPlot

The algorithm used by Plot samples functions at machine precision numbers.  In many cases this allows faster evaluation and almost always works very well.  Sometimes an acceptable graph can't be made using machine precision arithmetic.  When this is the case a non-trivial trick or two is needed to get Plot to produce an adequate graphic.  This is often the case if you want to make a plot to show how rapidly the difference between a function and it's asymptotic expansion approaches zero.  Consider the function (f1) below which is the difference between Exp[x^2] Erfc[x] and the (9^th) order asymptotic expansion of this function.

\(ClearAll[f1, x]; 

    f1[x_] = Exp[x^2] Erfc[x] - 
		
        Exp[x^2]*\((Series[Erfc[x], {x, ∞, 9}] // Normal)\)\)

If you evaluate the next cell on a Mac G3 using Version 3.0, you will get an erratic plot.  If you evaluate it with a Pentium using Version 3.0 you will get a smooth graph but the function values will have the wrong sign and the magnitude will be off by several orders of magnitude.  The results vary among different computers, but few if any computers will produce a satisfactory graphic using the line in the next cell.

Plot[f1[x], {x, 20, 30}];

The next cell shows how one can get the Plot algorithm to make a satisfactory graph of (f1) over (20<x<30).

ClearAll[g, h]g[x_] := f1[x]h[x_] := N[f1[SetPrecision[x, ∞]], 17]Plot[h[x], {x, 20, 30}];

The PrecisionPlot package defines a new function called PrecisionPlot to automate the code above.  Before you try the examples below you should put the PrecisionPlot.m package in the 3.0\AddOns\ExtraPackages folder and evaluate the next cell to load the package.

<< PrecisionPlot.m

The PrecisionPlot usage message is shown below.

?PrecisionPlot

Examples

In the next cell PrecisionPlot is used to make the same graphic that was made using Plot above.  When PrecisionPlot is used all the details of converting to arbitrary precision are hidden.  The thought process involved in writing the PrecisionPlot package is explained below in "The How and Why Behind PrecisionPlot".

ClearAll[f1]f1[x_] := 
      -(105/(16 (Sqrt[π]) (x^9))) + (15/(8 (Sqrt[π]) (x^7))) - 
        (3/(4 (Sqrt[π]) (x^5))) + (1/(2 (Sqrt[π]) (x^3))) - 
        (1/((Sqrt[π]) x)) + (E^((x^2))) Erfc[x]

PrecisionPlot[f1[x], {x, 20, 30}];

The cells that follow give another example where Plot doesn't work well.

\(ClearAll[x, f2]; 

    f2[x_] = ArcSec[x] - 
		
        Normal[Series[TrigToExp[ArcSec[x]], {x, ∞, 7}]]\)

Plot[f2[t], {t, 25, 40}];

PrecisionPlot uses the same options as Plot, and a RationalValues option that's unique to this package.  The usage message for this new option is shown below, and it has the default setting RationalValues->True.

?RationalValues

In the next cell PrecisionPlot produces a smooth graph when Plot couldn't in the last example.  Here the option RationalValues->False is used.

PrecisionPlot[f2[t], {t, 25, 45},  RationalValues -> False];

Other examples where Plot has trouble are given in the next three input cells.

\(ClearAll[x, f3]; 

    f3[x_] = Exp[1/x] - Normal[Series[Exp[1/x], {x, ∞, 8}]]\)

Plot[f3[t], {t, 10, 20}];

Plot[f3[t], {t, 20, 30}];

Next PrecisionPlot is used to make smooth graphs when Plot couldn't in the two previous cells.

PrecisionPlot[f3[t], {t, 10, 20}];

PrecisionPlot[f3[t], {t, 20, 30}];

In the next example PrecisionPlot is given a list of functions to plot and a different style for each function.  In addition the functions are plotted even though (t) evaluates to 3.  Special care was taken in writing the package to ensure PrecisionPlot works like Plot in such an example.  In fact PrecisionPlot will take any combination of options Plot would take and use them accordingly.

\(t = 3; 

    PrecisionPlot[{f3[t], 10\^\(-13\)\/\(4  t\), 10\^\(-14\)\/t}, 

      		{t, 10, 20}, Frame -> True, 
	\ \ 
      PlotStyle -> {{Hue[0.4], Thickness[0.0125]}, {Hue[0.7]}, {Hue[1]}}, 

      \ \ \ AspectRatio -> 1/2]; \)

The How and Why Behind PrecisionPlot

To be certain the cells that follow give the expected results you should evaluate the next cell.  The functions (f1) and (f3) are the same functions used above.

\(Unprotect[PrecisionPlot]; 
Remove[PrecisionPlot]\)Remove["Graphics`PrecisionPlot`Private`*"]\(ClearAll[f1, f3, t]
\ \)\(f1[x_] = 
	
      \(-\(105\/\(16\ \@π\ x\^9\)\)\) + 15\/\(8\ \@π\ x\^7\) - 
        3\/\(4\ \@π\ x\^5\) + 1\/\(2\ \@π\ x\^3\) - 
        1\/\(\@π\ x\) + E\^\(x\^2\)\ Erfc[x]; 


    f3[x_] = 
	
      \(-1\) + E\^\(1\/x\) - 1\/\(40320\ x\^8\) - 1\/\(5040\ x\^7\) - 
        1\/\(720\ x\^6\) - 1\/\(120\ x\^5\) - 1\/\(24\ x\^4\) - 
        1\/\(6\ x\^3\) - 1\/\(2\ x\^2\) - 1\/x; \)

We see why the Plot algorithm sometimes fails

Using the line in the next cell the Plot algorithm makes a smooth graph.

Plot[f3[t], {t, 4, 7}];

But the Plot algorithm makes an erratic graph using the line in the next cell.

Plot[f3[t], {t, 10, 20}];

The Plot algorithm samples a function at machine precision numbers.  Below we see that the function value is significantly different when it's evaluated at the machine precision number (18.0) compared to evaluating at (18) using arbitrary precision arithmetic.  If you want to ensure you results are correct you have to use arbitrary precision arithmetic.

f3[18.0]

N[f3[18], 17]

We look at methods of forcing Plot to use arbitrary precision arithmetic

Near the top of this notebook a method was given that will ensure Plot produces the right graphic.  Using that method the function given to Plot was evaluated at numbers with infinite precision.  Before the author discovered a way to quickly convert a real number to a rational number, two other methods were provided by Dave Withoff (of Wolfram Research) through technical support.  Each of these methods evaluate the function at arbitrary precision numbers having 17 digits of precision.  The first method from Dave Withoff is used in the next cell.  For reasons the author doesn't fully understand the option (Compiled->False) must be used with this method.

Plot[f3[SetPrecision[t,17]],{t,10,20},
Compiled->False,PlotRange->All];

The second method Dave Withoff provided requires the use of an intermediate function (h), as shown in the cells below.  When done this way the option (Compiled->True) or (Compiled->False) can be used.  This second method is used in the PrecisionPlot package.

ClearAll[h];
h[x_]:=f3[SetPrecision[x,17]];

Plot[h[t],{t,10,20},Compiled->True];

Plot[h[t],{t,10,20},Compiled->False];

We make a first attempt at writing a PrecisionPlot function

In the next cell we try to make a function PrecisionPlot by defining an intermediate function (h) in the cell below.  However, you can try it and see that it doesn't work.

\(PrecisionPlot[f_, {x_, xmin_, xmax_}, opts___?OptionQ] := \ 
      Module[{h}, 
\ \ \ \ 
        \((h = \((f[SetPrecision[#, 17]]&)\); 
\ \ \ \ \ \ \ 
          Plot[h[x], {x, xmin, xmax}, opts]
		)\)]\)

The attempt above fails because the Plot algorithm tries to sample the function (f[SetPrecision[#,17]]&) but the kernel doesn't know (x) is the argument of (f).  This is demonstrated in the cell below.

\(ClearAll[f, x]; 
f = x + Exp[x]; 
f[5]\)

Instead we can make a pure function (g) where each (x) in (f) becomes the argument of (g) as in the next cell.

\(ClearAll[g, x]; 
g = Function[f /. x -> #]; 
g[5]\)

We make our first working PrecisionPlot function

Our first working version of PrecisionPlot is given below.  We use the form Function[f/.x->#] instead of the equivalent form f/.x->#& to ensure the definition is more easily understood after some modifications are made later in this notebook.

\(ClearAll[PrecisionPlot]; 

    PrecisionPlot[f_, {x_, xmin_, xmax_}, opts___?OptionQ] := \ 
      Module[{g, h}, 
\ \ \ \ 
        \((g = Function[f /. x -> #]; 
\ \ \ \ \ \ \ 
          h = \((g[SetPrecision[#, 17]]&)\); 
\ \ \ \ \ \ \ 
          Plot[h[x], {x, xmin, xmax}, opts]
		)\)]\)

In the next cell PrecisionPlot is used to make a plot that requires the use of arbitrary precision arithmetic.

ClearAll[t];
PrecisionPlot[f3[t],{t,10,20}];

However, PrecisionPlot works poorly in the next example where the option (Compiled->False) is used.

PrecisionPlot[f3[t], {t, 10, 20}, Compiled -> False];

We ensure PrecisionPlot works when using the setting (Compiled->False)

The problem with the last definition of PrecisionPlot can be corrected by using (g=Function[Evaluate[f/.x->#]]) instead of (g=Function[f/.x->#]) in the definition of PrecisionPlot.  This solution was also provided by David Withoff via technical support and is used in the next cell.  However, the author is unable to explain why the previous definition for PrecisionPlot works well using (Compiled->True) but not using (Compiled->False).

\(ClearAll[PrecisionPlot]; 

    PrecisionPlot[f_, {x_, xmin_, xmax_}, opts___?OptionQ] := \ 
      Module[{g, h}, 
\ \ \ \ 
        \((g = Function[Evaluate[f /. x -> #]]; 
\ \ \ \ \ \ \ 
          h = \((g[SetPrecision[#, 17]]&)\); 
\ \ \ \ \ \ \ 
          Plot[h[x], {x, xmin, xmax}, opts]
		)\)]\)

Using this new definition for PrecisionPlot the next two cells produce satisfactory graphics using (Compiled->False) and (Compiled->True).

PrecisionPlot[f3[t], {t, 10, 20}, Compiled -> False];

PrecisionPlot[f3[t], {t, 10, 20}, Compiled -> True];

We give PrecisionPlot Options and Attributes

Then we add a few extras in the next cell.
  (1) We give PrecisionPlot the same options as Plot.  
       We also give it the default setting Compiled->False since CompiledFunctions 
       can't use arbitrary precision arithmetic.  However, using the setting Compiled->True 
       makes no difference of practical significance.
  
  (2) We ensure the definition isn't used to plot a list of functions.
  
  (3) We give PrecisionPlot the HoldAll attribute so it will work the same as 
         Plot in this regard.  In the package PrecisionPlot is also protected.
         
Evidence isn't provided that the three features above are correctly implemented, but you can easily do this if you like.

\(ClearAll[PrecisionPlot]; 

SetAttributes[PrecisionPlot, HoldAll]\)\(Options[PrecisionPlot] = Options[Plot]; 

    SetOptions[PrecisionPlot, Compiled -> False]; 


    PrecisionPlot[f_, {x_, xmin_, xmax_}, 
opts___?OptionQ] /; 
        \((Head[f] =!= List)\) := \ 
\ 	
      Module[{g, h}, 
\ \ \ \ \ g = Function[Evaluate[f /. x -> #]]; 

        \ \ \ \ \ h = g[SetPrecision[#, 17]]&; 
\ \ \ \ \ 
        Plot[h[x], {x, xmin, xmax}, opts]
\ \ ]; \)

We ensure a local environment is used for the dummy variable

The PrecisionPlot definition developed in the previous group of cells is copied below.  If this isn't the current definition of PrecisionPlot you should evaluate the next cell at this time.

\(ClearAll[PrecisionPlot]; 

SetAttributes[PrecisionPlot, HoldAll]\)\(Options[PrecisionPlot] = Options[Plot]; 

    SetOptions[PrecisionPlot, Compiled -> False]; 


    PrecisionPlot[f_, {x_, xmin_, xmax_}, 
opts___?OptionQ] /; 
        \((Head[f] =!= List)\) := \ 
\ 	
      Module[{g, h}, 
\ \ \ \ \ g = Function[Evaluate[f /. x -> #]]; 

        \ \ \ \ \ h = g[SetPrecision[#, 17]]&; 
\ \ \ \ \ 
        Plot[h[x], {x, xmin, xmax}, opts]
\ \ ]; \)

Consider the next cell where the PrecisionPlot algorithm above is used to make a graphic that doesn't require the use of arbitrary precision.

\(ClearAll[t]; 
PrecisionPlot[\@12 + Log[t] + t\^2/100, {t, 1, 80}]; 
    \)

In the next cell Plot makes the same graph even though (t) has a global value.

\(t = 3; 
Plot[\@12 + Log[t] + t\^2/100, {t, 1, 80}]; \)

Then we try to do the same thing using PrecisionPlot, and we get a graph that isn't even close to the expected graph.  This is very undesirable since one could unknowingly make a graph of the wrong function.

\(t = 3; 
PrecisionPlot[\@12 + Log[t] + t\^2/100, {t, 1, 80}]; \)

The problem demonstrated in the previous cell is due to expression (g=Function[Evaluate[f/.x->#]]) in the definition of PrecisionPlot.  Through technical support Dale Horton (of Wolfram Research) reported that this problem can be solved using Block or Module.  To demonstrate consider the function (fun) below.

\(ClearAll[fun]; 
SetAttributes[fun, HoldAll]; 

    fun[expr_, sym_] := Function[Evaluate[expr /. sym -> #]]\)

In the cell below (t) has no values and each occurrence of (t) is changed to (#1) as intended.

\(ClearAll[t]; 
fun[\@12 + Log[t] + t^2/100, t]\)

Now in the next line (t=3) and we don't get the intended result.

\(t = 3; 
fun[\@12 + Log[t] + t^2/100, t]\)

The expression in the last cell didn't give the intended result because the evaluation process went as shown below. Wrapping the expression in Trace will reveal this.
 
fun[(Sqrt[12]) + Log[t] + (t^2)/100, t]
=> Function[Evaluate[(Sqrt[12]) + Log[t] + (t^2)/100  /. t -> #]]
=>(2 (Sqrt[ 3]) + Log[t] + (t^2)/100  /. t -> #1)&
=>(2 (Sqrt[ 3]) + Log[3] + (t^2)/100  /. t -> #1)&
=> (2 (Sqrt[ 3]) + Log[3] + 9/100  /.  t -> #1)&
=> (9/100 + 2 (Sqrt[ 3]) + Log[3]  /.  t -> #1)&
=> (9/100 + 2 (Sqrt[ 3]) + Log[3]  /.  3 -> #1)&
=> \(\((9/100 + 2 \@ #1 + Log[#1])\)&\)
=> \(\((9/100 + Log[#1] + 2 \@ #1)\)&\)
The important point in the evaluation above is that (t) evaluates to (3) and ((Sqrt[12])) evaluates to (2 (Sqrt[ 3])) before the replacement is made.  When the replacement is made the rule (3->#1) is used.  A different function (funModule) is given in the next cell where (sym) is put in a module.

\(ClearAll[funModule]; 
SetAttributes[funModule, HoldAll]; 

    funModule[expr_, sym_] := 
      Function[Evaluate[Module[{sym}, expr /. sym -> #]]]\)

In the next cell funModule replaces each occurrence of (x) with (#1) even though (x=3).

\(t = 3; 
funModule[\@12 + Log[t] + t^2/100, t]\)

When funModule evaluated in the last example the kernel selected a unique symbol (e.g.  t$58) to use instead of (t).  Supposing the unique symbol was (t$58) the evaluation process for funModule above went as shown below.  Wrapping the expression in Trace will reveal this as well.

funModule[(Sqrt[12]) + Log[t] + (t^2)/100, t]
=> Function[Evaluate[Module[{t}, (Sqrt[12]) + Log[t] + (t^2)/100  /. t -> #]]]
        
=>((Sqrt[12]) + Log[t$58] + (t$58^2)/100  /. t$58 -> #1)&
=>(2 (Sqrt[ 3]) + Log[t$58] + (((t$58))^2)/100  /. t$58 -> #1)&
=> (2 (Sqrt[ 3]) + (((t$58))^2)/100 + Log[t$58]  /.  t$58 -> #1)&
=> (2 (Sqrt[ 3]) + (((#1))^2)/100 + Log[#1])&
=> (2 (Sqrt[ 3]) + Log[#1] + (((#1))^2)/100)&
This time every occurrence of (t) was changed to the unique symbol before the subexpressions evaluated and the replacement was done as intended.  Alternatively we can use Block instead of Module to ensure the replacement is done correctly.  The function funBlock is defined below to demonstrate this.

\(ClearAll[funBlock]; 
SetAttributes[funBlock, HoldAll]; 

    funBlock[expr_, sym_] := 
      Function[Evaluate[Block[{sym}, expr /. sym -> #]]]\)

In the next cell funBlock replaces each occurrence of (t) with (#1) even though (t=3).

\(t = 3; 
funBlock[\@12 + Log[t] + t^2/100, t]\)

The evaluation process for funBlock in the last example is shown below.
 
\(TextForm
      \`funBlock[\@12 + Log[t] + t\^2/100, t]
 => 
        Function[Evaluate[Block[{t}, \@12 + Log[t] + t\^2/100\  /. t -> #]]]
          
 => 
        Values\ of\ \((t)\)\ are\ temporarily\ \(cleared . \)\)
 => (2 (Sqrt[ 3]) + Log[t] + (t^2)/100  /.  t -> #)&
 => (2 (Sqrt[ 3]) + (t^2)/100 + Log[t]  /. t -> #)&
 => 2 (Sqrt[3]) + #((1^2)/100) + Log[#1]&
 =>2 (Sqrt[3]) + Log[#1] + #((1^2)/100)&
 => Original values of (t) are restored.
 
In this case the global values of (t) were cleared before the subexpressions were evaluated, and the replacement was made as intended.  After making the replacement Block restored the original global values to (t).

When funModule evaluates the kernel has to replace each occurrence of (t) with the unique symbol.  There may be many occurrences of (t) in the given expression, and funBlock would be considerably faster in such a case since it doesn't have this extra overhead.  For this reason we use Block to correct the last problem with PrecisionPlot.  The new definition of PrecisionPlot is shown below.

\(ClearAll[PrecisionPlot]; 

SetAttributes[PrecisionPlot, HoldAll]\)\(Options[PrecisionPlot] = Options[Plot]; 

    SetOptions[PrecisionPlot, Compiled -> False]; 


    PrecisionPlot[f_, {x_, xmin_, xmax_}, opts___?OptionQ] /; 
	
        \((Head[f] =!= List)\) := \ 
\ 	
      Module[{g, h}, 
		\ \ 
        g = Function[Evaluate[Block[{x}, f /. x -> #]]]; 
\ \ \ \ \ \ \ 
        h = g[SetPrecision[#, 17]]&; 
\ \ \ \ \ \ \ 
        Plot[h[x], {x, xmin, xmax}, opts]
\ \ ]; \)

The new version of PrecisionPlot is demonstrated below and makes the intended graph when (t) has a numeric value.

\(t = 3; 
PrecisionPlot[\@12 + Log[t] + t^2/100, {t, 1, 81}]; \)

We force PrecisionPlot to evaluate functions at rational numbers

The PrecisionPlot definition developed in the previous group of cells is copied below.  If this isn't the current definition of PrecisionPlot you should evaluate the next cell at this time.

\(ClearAll[PrecisionPlot]; 

SetAttributes[PrecisionPlot, HoldAll]\)\(Options[PrecisionPlot] = Options[Plot]; 

    SetOptions[PrecisionPlot, Compiled -> False]; 


    PrecisionPlot[f_, {x_, xmin_, xmax_}, opts___?OptionQ] /; 
	
        \((Head[f] =!= List)\) := \ 
\ 	
      Module[{g, h}, 
		\ \ 
        g = Function[Evaluate[Block[{x}, f /. x -> #]]]; 
\ \ \ \ \ \ \ 
        h = g[SetPrecision[#, 17]]&; 
\ \ \ \ \ \ \ 
        Plot[h[x], {x, xmin, xmax}, opts]
\ \ ]; \)

There is yet another problem with the program in the cell above.  Andrzej Kozlowski pointed out that the program developed above doesn't work well in the next example.

PrecisionPlot[f1[x], {x, 70, 80}];

Using Version 3.0 on a Mac G3 or a PC the cell above produces a graphic with a discontinuity at about (75.2). The function f1[x] (defined earlier) is clearly continuous, so the PrecisionPlot algorthm doesn't do well in this case.  To see what went wrong we consider what the function value is at an arbitrary precision number having 17 digits of precision.

\(tst = SetPrecision[76.0, \ 17]; 
y = f1[tst]\)

The output of the previous cell does nothing to help explain the graph in the previous example.  To understand what happened with the previous graph consider the InputForm of the result above.

\($NumberMarks = True; 
InputForm[y]\)

The InputForm above is much more helpful.  We see the function evaluates to "approximately" -1.3495009×(10^(-19)) with -4.0469 digits of precision.  A better estimate of what (f1) evaluates to is given in the cell below.  This is considerably different from the first value.  Given that f1[tst] had negative precision we shouldn't be surprised that the result wasn't very close to the correct value.  The incorrect values were used to make the plot with no condideration to the fact that they had negative precision.

N[f1[76],  17]

In order to ensure that PrecisionPlot will compute the funtion values with sufficient precision we need to use arbitrary precision arithmetic to evaluate the function at numbers with infinite precision.  The Plot algorithm samples functions at values with the head Real.  We could use Rationalize to convert to rationals, but this is relatively slow.  SetPrecision has an undocumented capability to convert to rationals as demonstrated below.

a = SetPrecision[74.8,  ∞]

The result above is obviously not as consise as that you would get from Rationalize, but SetPrecision (as used above) evaluates about 26 times fater than Rationalize[74.8, 0.].  The definition below will ensure PrecisionPlot evaluates functions at rational numbers with sufficient precision.

\(ClearAll[PrecisionPlot]
\)SetAttributes[PrecisionPlot, HoldAll]\(Options[PrecisionPlot] = Options[Plot]; 

    SetOptions[PrecisionPlot, Compiled -> False]; 


    PrecisionPlot[f_, {x_, xmin_, xmax_}, opts___?OptionQ] /; 
	
        \((Head[f] =!= List)\) := \ 
\ \ 
      Module[{g, h}, 
\ \ \ \ 
        g = Function[Evaluate[Block[{x}, f /. x -> #]]]; 
		
        h = N[g[SetPrecision[#, Infinity]], 17]&; 
\ \ \ \ 
        Plot[h[x], {x, xmin, xmax}, opts]
\ ]; \)

For many purposes the relatively short program above is all you need.  Similar programs can be written that use arbitrary precision for Plot3D, ParametricPlot, ContourPlot, etc.  In the line below we demonstrate that this new definition for PrecisionPlot corrects the problem demonstrated earlier.  In this case PrecisionPlot produces an aceptable plot of f1[x].

PrecisionPlot[f1[x], {x, 70, 80}];

We give PrecisionPlot a RationalValues Option

Some users might have a problem with the definition of PrecisionPlot developed in the previous group of cells.  This latest definition often takes about three times longer than the ealier definition of PrecisionPlot.  To address this concern we give PrecisionPlot an option that allows the user to control which method is used.

This is done using the code in the next cell.  Notecie the FilterOptions package is used to ensure only the appropriate options are passed to Plot.  In addition the options passed to Plot must be evaluated before the Plot algorithm begins.

Needs["Utilities`FilterOptions`"]\(ClearAll[PrecisionPlot]
\)\(SetAttributes[PrecisionPlot, HoldAll]
\)\(Options[PrecisionPlot] = 
	
      Join[{RationalValues -> True}, Options[Plot]]; 


    SetOptions[PrecisionPlot, Compiled -> False]; 


    PrecisionPlot[f_, {x_, xmin_, xmax_}, opts___?OptionQ] /; 
\ \ 
        \((Head[f] =!= List)\) := \ 
\ \ 
      Module[{g, h, rv, opts1}, 
		
        rv = RationalValues /. Flatten[{opts, Options[PrecisionPlot]}]; 
		
        opts1 = FilterOptions[Plot, opts]; 
		
        g = Function[Evaluate[Block[{x}, f /. x -> #]]]; 
		
        If[rv, 
		\ \ h = N[g[SetPrecision[#, Infinity]], 17]&, 
		\ \ 
          h = g[SetPrecision[#, 17]]&
	\ ]; 
\ \ \ \ 
        Plot[h[x], {x, xmin, xmax}, Evaluate[opts1]]
\ ]; \)

Using the the next two cells PrecisionPlot makes essentialy the same graph.  When the option (RationalValues->False) is used the graph in made in about one third the time it takes using (RationalValues->True).

PrecisionPlot[f1[x], {x, 40, 50}];

PrecisionPlot[f1[x], {x, 40, 50}, RationalValues -> False];

However (using Version 3.0 on a PC) only the first of the next two lines produce the right graph. Also when the option (RationalValues->False) is used the program takes longer to make the graph because it samples the function many times where the function seems erratic.  The bottom line is that you risk getting a misleading graph if you use the setting (RationalValues->False).

PrecisionPlot[f1[x], {x, 76, 86}];

PrecisionPlot[f1[x], {x, 76, 86}, RationalValues -> False];

As another example consider the next cell where (using Version 3.0 on a Personal Computer) PrecisionPlot produces a smooth graph in about a third the time it takes using the option (RationalValues->True), but unbeknownst to the user the function values are in error by a noticable amount.

PrecisionPlot[f1[x], {x, 76, 106}, RationalValues -> False];

PrecisionPlot[f1[x], {x, 76, 106}];

We ensure PrecisionPlot can take a list of functions with a single PlotStyle

As indicated earlier the PrecisionPlot function developed thus far can't be used to plot a list of functions.  However, we can make provisions so that when PrecisionPlot is given a list of functions {f1,f2,f3, ...} the above PrecisionPlot algorithm is mapped to each (fi).  This approach will work as long as the PlotStyle option doesn't specify a different style for different functions.  A function called MultiplePlotStyles is defined below to address this possibility, and returns True if it's arguments imply a list of lists for the PlotStyle option.

\(ClearAll[MultiplePlotStyles]; 


    MultiplePlotStyles[plotFunc_, opts___] := \ 
	
      With[{
	\ \ \ 
          ps = PlotStyle /. 
	\ \ \ Flatten[{opts, Options[plotFunc]}]

          \ \ }, 
\ \ MatchQ[ps, {_List, __List}]
]; \)

If we want to evaluate 
PrecisionPlot[{f1,f2,f3},{x,xmin,xmax}]
where the default DisplayFunction is $DisplayFunction we need a PrecisionPlot definition that will use the simpler version of PrecisionPlot as follows:

Show[{
PrecisionPlot[f1,{x,xmin,xmax},DisplayFunction->Identity],
PrecisionPlot[f2,{x,xmin,xmax},DisplayFunction->Identity],
PrecisionPlot[f3,{x,xmin,xmax},DisplayFunction->Identity]
},DisplayFunction->$DisplayFunction]

Here PrecisionPlot is mapped to each (fi) using (DisplayFunction->Identity) so the graphics are not rendered.  Then all the graphics are sent to Show using (DisplayFunction->$DisplayFunction) and the graphics are rendered.  The PrecisionPlot definition below does this with the condition Not[MultiplePlotStyles[PrecisionPlot,opts]] to ensure the PlotStyle is not a list of lists.

Another detail in the definition below is the use of  Block[{x}, .....].  This is used because (f) and (x) evaluate before the simpler version of PrecisionPlot is used, and we want to ensure global values of the symbol represented as (x) have no effect on the graphic produced.

\(PrecisionPlot[f : {__}, {x_, xmin_, xmax_}, 
	opts___?OptionQ] /; 

        \ \ \ \ \ \ Not[MultiplePlotStyles[PrecisionPlot, opts]] := \ 

      With[
	{
\ \ \ 	
          df = DisplayFunction /. 
	\ \ \ \ 
              Flatten[{opts, Options[PrecisionPlot]}]
\ \ }, 
\ \ 
        Show[
		\ 
          Block[{x}, 
		\ \ \ \ \ 
            Map[
				\ \ \ 
              PrecisionPlot[#, {x, xmin, xmax}, 
\ \ \ \ \ \ \ \ \ \ \ \ 
                  DisplayFunction -> Identity, opts]&, f
			\ \ ]

            	\ \ \ \ ], DisplayFunction -> df
\ \ \ ]
]\)

Using the last definition for PrecisionPlot we get satisfactory graphics for each of the following examples.  Since Block[{x}, ...] was used in the last definition for PrecisionPlot the graphic produced isn't effected by global values of (t).

\(t = 3; 

    PrecisionPlot[{f3[t], 10\^\(-13\)\/\(4  t\), 10\^\(-14\)\/t}, {t, 10, 
        30}]; \)

\(\(PrecisionPlot[{f3[t], 10\^\(-13\)\/\(4  t\), 10\^\(-14\)\/t}, {t, 10, 
        20}, 
PlotStyle -> Automatic]; \)\)

\(\(PrecisionPlot[{f3[t], 10\^\(-13\)\/\(4  t\), 10\^\(-14\)\/t}, {t, 10, 
        20}, 
PlotStyle -> {}]; \)\)

\(\(PrecisionPlot[{f3[t], 10\^\(-13\)\/\(4  t\), 10\^\(-14\)\/t}, {t, 10, 
        20}, 
PlotStyle -> Dashing[{0.02, 0.02}]]; \)\)

\(\(PrecisionPlot[{f3[t], 10\^\(-13\)\/\(4  t\), 10\^\(-14\)\/t}, {t, 10, 
        20}, 
PlotStyle -> {Dashing[{0.02, 0.02}], Hue[0.7]}]; \)\)

\(\(PrecisionPlot[{f3[t], 10\^\(-13\)\/\(4  t\), 10\^\(-14\)\/t}, {t, 10, 
        20}, 
PlotStyle -> {{Dashing[{0.02, 0.02}], Hue[0.7]}}]; \)\)

\(\(PrecisionPlot[{f3[t], 10\^\(-13\)\/\(4  t\), 10\^\(-14\)\/t}, {t, 10, 
        20}, 
	Axes -> False, {
        PlotStyle -> Dashing[{0.02, 0.02}]}, 
AspectRatio -> 1/2]; 
    \)\)

\(\(PrecisionPlot[{f3[t], 10\^\(-13\)\/\(4  t\), 10\^\(-14\)\/t}, {t, 10, 
        20}, 
	Axes -> False, {
        PlotStyle -> {Dashing[{0.02, 0.02}], Hue[0.7]}}, 

      AspectRatio -> 1/2]; \)\)

We ensure PrecisionPlot can take a list of functions with multiple PlotStyles

We still need to address the case where PrecisionPlot is used to Plot a list of functions and a different style is used for different curves.  The next two examples demonstrate that when using Plot the number of styles given by PlotStyle need not be equal to the number of functions plotted.  In the cells below the Plot algorithm uses one of the styles for each function starting with the first style.  If necessary, the algorithm cycles through the list of styles again and again.

\(\(Plot[{Sin[t], Sin[2  t], Sin[3  t], Sin[4  t], Sin[5  t]}, 
      
{t, 0, π/6}, PlotStyle -> {{Hue[1]}, {Hue[0.3]}}]; \)\)

\(\(Plot[{Sin[4  t], Sin[5  t]}, {t, 0, π/6}, 
      PlotStyle -> 

        		{{Hue[1]}, {Hue[0.3]}, {Hue[0.6]}, {Hue[0.9]}}]; \)\)

We want to ensure PrecisionPlot can do the same.  Towards this end we define a function BuildPlotStyles which takes a list of functions and a list of styles and returns a list of styles with the same length as the list of functions.

\(ClearAll[BuildPlotStyles]; 


    BuildPlotStyles[functs_List, styles_List] := 

      Take[Join@@
          Table[styles, 
\ \ \ \ 	{Ceiling[Length[functs]/Length[styles]]}

            	\ ], 
\ 	{1, Length[functs]}
]\)

The next two cells provide demonstration of the BuildPlotStyles function.

\(BuildPlotStyles[{x, 2  x, 3  x, 4  x, 5  x, 6  x, 7  x}, 

      	{{opts1 -> True}, {opts2 -> True, opts3 -> False}, {
          opts4 -> True}}]\)

\(BuildPlotStyles[{x, 2  x, 3  x}, 

      	{{opts1 → True}, {opts2 → False, 
          opts3 -> True}, {opts4 → True}, 

        		{opts5 → True}, {opts6 → True}, {
          opts7 → True}, 
		{opts8 → True}}]\)

Now if we want to evaluate 
PrecisionPlot[{f1,f2,f3},{x,xmin,xmax}, PlotStyle->{{style1},{style2},{style3}}]
where the DisplayFunction used is $DisplayFunction we can ensure that the following expression is used to achieve this.

Show[{
  PrecisionPlot[f1,{x,xmin,xmax}, DisplayFunction->Identity, PlotStyle->style1],
  PrecisionPlot[f2,{x,xmin,xmax}, DisplayFunction->Identity, PlotStyle->style2],
  PrecisionPlot[f3,{x,xmin,xmax}, DisplayFunction->Identity, PlotStyle->style3]
},DisplayFunction->$DisplayFunction]

Here MapThread can be used to evaluate the simpler version of PrecisionPlot with the respective functions and styles.  As with the previous definition the simpler version of PrecisionPlot is evaluated with the option (DisplayFunction->Identity) to ensure graphics are not rendered.  Then all the graphics are sent to Show using (DisplayFunction->$DisplayFunction) and the graphs are rendered.

Here again (f) and (x) are evaluated inside Block[{x}, ....] to ensure global values of the symbol represented as (x) have no effect on the graphic produced.

\(PrecisionPlot[f : {__}, {x_, xmin_, xmax_}, opts___?OptionQ] /; 
	
        MultiplePlotStyles[PrecisionPlot, opts] := 

      Module[{df, ps, graph}, 

        \ \ {ps, df} = {PlotStyle, DisplayFunction} /. 
	\ \ \ 
            Flatten[{opts, Options[PrecisionPlot]}]; 
\ \ \ \ \ 
        Block[{x}, \ 
		\ \ \ \ ps = BuildPlotStyles[f, ps]; 

          \ \ \ \ \ \ \ \ 
          graph = MapThread[
			\ \ \ \ \ 
              PrecisionPlot[
		\ \ \ \ \ \ \ \ \ \ #1, {x, xmin, xmax}, 
                  DisplayFunction -> Identity, 
			\ \ \ \ \ \ \ \ 
                  PlotStyle -> #2, opts
			\ \ \ \ ]&, {f, ps}

              	\ \ \ \ \ ]; 
	]; 
\ \ \ 
        Show[graph, DisplayFunction -> df]
]\)

Three examples where PrecisionPlot uses a different style for different curves are given in the cells below, and the graphics are not effected by global values of (t).

\(t = 3; 

    PrecisionPlot[{f3[t], 10\^\(-13\)\/\(4  t\), 10\^\(-14\)\/t}, {t, 10, 
        20}, 
PlotStyle 
        -> {{Hue[0.3], Thickness[0.02]}, {Hue[0.7]}, {Hue[1]}, {
            Dashing[{0.02, 0.02}]}}]; \)

\(\(PrecisionPlot[{f3[t], 10\^\(-14\)\/t}, {t, 10, 20}, 

      PlotStyle 
        -> {{Hue[0.3], Thickness[0.02]}, {Hue[0.7]}, {Hue[1]}, {
            Dashing[{0.02, 0.02}]}}]; \)\)

\(\(PrecisionPlot[{f3[t], 10\^\(-13\)\/\(4  t\), 10\^\(-13\)\/\(8\ t\), 
        10\^\(-14\)\/t}, {t, 10, 20}, 

      PlotStyle -> {{Hue[0.3], Thickness[0.013]}, {Hue[1]}}]; \)\)

We put the definitions in a package

Finally the above definitions are put in a package using the usual package constructs.  In particular the code below is used to achieve this.

BeginPackage["Graphics`PrecisionPlot`", "Utilities`FilterOptions`"]
PrecisionPlot::usage="PrecisionPlot[f, {x, xmin, xmax}]....."
Begin["`Private`"]
   .....
   ..... (* PrecisionPlot definitions *)
   ..... (* MultiplePlotStyles definition *)
   ..... (* BuildPlotStyles definition *)
   .....
   Protect[PrecisionPlot];
End[]
EndPackage[];