This manual describes MICROSCOPE: a portable FORTRAN software system for the analysis of multivariate functions. Given an interpolation or approximation scheme, it allows the following questions, among others, to be answered:
This manual constitutes the comprehensive description of MICROSCOPE, its operation, applications, foundations, and inner workings.
Section 1.1 describes what kind of problems can be addressed by MICROSCOPE. Section 1.2 describes the basic approach and is fundamental. Section 1.3 describes the alphanumerical screen display and is needed for understanding the examples given throughout this manual. It also introduces some terms and concepts and serves as a prelude to the capabilities of MICROSCOPE.
Section 2 gives examples that illustrate four major applications (interpolation, smoothness, degree of precision, and polynomial degree). These do of course not exhaust the potential of MICROSCOPE. In section 3, one other application (discovery or identification of a point where smoothness is reduced) is used as a vehicle for introducing the mechanics of using MICROSCOPE.
Section 4 gives a detailed description of all commands and features, organized by type of feature or task to be accomplished. In regular use, this will probably be the most frequently consulted section.
In reading sections 2 and 3, it is desirable that a working version of MICROSCOPE is available so that hands on experience can be gained. If this is not the case the installation guide (section 5) should be consulted first. That section also contains a detailed discussion of the portability features and possible variants of MICROSCOPE.
Appendices are included that describe some technical features in detail: the differentiation formulas, the organization of the program package, and the test package. There should also be an appendix IV containing pen plots corresponding to the cruder alphanumerical graphics in the body of the manual. However, that appendix is not machine legible, and is distributed separately.
The manual itself is in machine readable form so that it can be distributed together with the package and can be printed locally, and read on line using for example the search features of a text editor. It was generated by a MICROSCOPE program, with the text of the manual being "notes" and the example screen displays being computed by MICROSCOPE and automatically embedded into the text. (Actually, the MICROSCOPE program generated a data file that was then fed into the text formatting program DOCUMENT, see Beebe, 1980.) Machine legibility, however, accounts for the somewhat awkward notation in some places. For example, greek letters have been given abbreviations in terms of roman letters.
Definitions are denoted by writing the defined term in single quotation marks. Ordinary quotations are in double quotation marks.
The examples in this manual are reproducible since the set of routines used for their generation is supplied with the package and can be used for test purposes. To ensure uniformity of results obtained on different machines, a facility has been added to simulate rounding, and all examples have been computed in 10 digits arithmetic. More precisely, if x is the quantity to be rounded to D digits, say, then x is replaced by -D -D z := (1+eps *10 )*x + eps *10 1 2 where eps and eps are random numbers between -1 and +1. The addition of 1 2 the eps term is not standard but appropriate in the present context, 2 because in investigations with MICROSCOPE small numbers are often due to taking differences between very close large numbers, leading to a cancellation of significant digits. The larger the magnitude of x, the more insignificant the second term will become. If x is roughly 1 or larger (in magnitude) our rounding is essentially equivalent to the standard rounding (where eps = 0). 2
This manual will rarely be read sequentially from front to back. Therefore, key ideas and concepts have sometimes been restated if they are crucial in the given context. It is hoped that this redundancy will contribute to the usefulness of this document.
Any critcisms, comments and suggestions about MICROSCOPE are welcome and should be directed to Peter Alfeld, Department of Mathematics, University of Utah, Salt Lake City, Utah 84112, 801-581-6842. We are particularly interested in any applications of MICROSCOPE that are not described in this manual.
This manual describes a FORTRAN software package, MICROSCOPE, whose main applications are in the analysis of functions of several variables (or 'surfaces') as they occur e.g. in the approximation of multivariate data or in the design of geometric objects such as the body of an automobile or an aircraft. This is a currently very active research area with many unsolved and difficult problems. For a survey see Barnhill, 1983, and the references quoted therein.
MICROSCOPE is portable to the extent that it has passed the PFORT verifier (Ryder, 1974) with the exception that it uses some non-FORTRAN characters (namely ?,<,>,!,:,;). Some non-portable optional features are supported by variants of the code. These include the use of lower case letters in Input and Output, and the use of modifications of an alphanumerical display on a CRT terminal. If the software ( <PLOT79>, Beebe 1979) and hardware are available more sophisticated and pleasing graphical displays can also be obtained. However, even if none of the additional features are provided by a particular installation, a working, if somewhat crude, version of MICROSCOPE can be constructed. See the installation guide (section 5) for details.
The objects that can be examined by MICROSCOPE are functions of from one to three independent variables. Allowing only up to three variables enabled us to contain all numerical information on the CRT screen and proved sufficent for all cases we have encountered. If a function of more than three variables must be examined MICROSCOPE can be applied to a suitable restriction of the function to a subdomain of dimension at most three. If there is sufficent interest, future versions of MICROSCOPE capable of handling more than three variables will be provided.
The functions of interest are usually obtained by applying an interpolation or approximation scheme to discrete or transfinite (i.e. infinitely many) data. As far as the use of MICROSCOPE is concerned it matters little if a function is obtained by interpolation or approximation (except that one will not need to test for interpolation in the latter case), so we will use the two terms interchangeably.
We will refer to the function that is being examined as the 'trial function'. The trial function may be different from the interpolating or approximating function in which the final user of a numerical scheme is interested. For example, it is often useful to generate data from a 'primitive function' and then use the difference between the primitive function and its interpolant as the trial function.
This makes discontinuities in derivatives large relative to the derivative value and thereby amplifies discontinuities which are then more easily detected.
Four types of questions occur most frequently in the analysis of the trial function:
It should be obvious that only in the very simplest of cases answers to these questions can be obtained from a plain display of the surface of interest. Even given infinite display accuracy (although of course any physical display is ultimately piecewise linear or even discrete), it seems impossible to distinguish optically a once differentiable function from one that is twice differentiable.
In the design of an approximation scheme one usually knows the 'critical sets' where the phenomena of interest such as discontinuities in derivatives, or interpolation to certain data, may occur. For example, in a bivariate interpolation scheme defined on a triangulation, critical sets are usually edges of triangles, possibly internal edges, and vertices of triangles. In a trivariate scheme defined on a tesselation into tetrahedra, critical sets are faces, edges, maybe internal faces and edges, and internal and external vertices. If the critical set is not known it can be discovered with MICROSCOPE (see section 3), but for the present we will assume it is known.
The basic idea of MICROSCOPE is very simple: Pick a random point (the 'point of examination') in a specific critical set. We will refer to that set as 'the front'. Then determine a 'direction of investigation' and consider a line through the point of examination in the direction of investigation. This is the 'line of investigation'. Depending upon the application, the line of investigation may be contained in the front, or it may be at an angle (not necessarily a right one) to it. We will call any derivative in the direction of investigation a 'tangential derivative'. Any pure derivative in a direction other than the direction of investigation is a 'cross derivative' and its direction is the 'cross direction'. MICROSCOPE numerically approximates the relevant derivatives (pure or mixed) and displays the approximations on the CRT terminal. Depending on the resulting picture, one may then take further action, e.g. pick another point of examination or another direction of investigation, consider other derivatives, or improve the existing plot by altering the underlying numerical parameters.
In using MICROSCOPE it is important to choose the underlying domain and the primitive function so as to avoid artifacts due to the regularity of the domain or peculiarities of the primitive function, and to pick the points of examination randomly or arbitrarily within the critical sets. Usually different types of critical sets will be present and of course one will need to consider at least one representative of each set.
Numerical differentiation is a notoriously ill-conditioned numerical process. However, in MICROSCOPE, the appearance of round-off errors is immediately apparent because the points on the graph of the relevant derivative are scattered across the screen in a random fashion. Moreover, by plotting derivatives of a high degree with a large discretization parameter the limitations imposed by round-off errors can usually be overcome. For example, if the expected discontinuity in a specific derivative is too unpronounced to be visible at the large value of the discretization parameter necessitated by round-off errors, then higher order derivatives can be computed (on an even coarser discretization). These will be approximations of the Dirac Delta function and its derivatives, and usually will pinpoint even weak discontinuities.
Numerically, we proceed as follows. Let F denote the trial function, P the point of examination, and d the direction of investigation. MICROSCOPE displays derivatives of F in the direction of d along the line of investigation with P being the center point. A built-in facility allows the replacement of F with a derivative of F in some other, independent, direction. More complicated mixed partial derivatives can be built by modifying the trial function, and MICROSCOPE offers some help with that as well. However, for the present we will restrict our attention to a tangential derivative of F.
The directional derivatives are approximated at points
P(i) = P + i*s*d, where i = -N, -N+1, ...., N-1,N
and the quantity s > 0 determines the spacing of the points. It will be refered to as the 'display interval' or 'interval' for short. The number N determines the resolution of the display (as well as the numerical effort). At each point P(i), we approximate the k-th derivative of F in the direction d at the i-th point by a quantity
(k) Q . i
The range of k in the present version of MICROSCOPE is from zero to six. We have found this sufficient for all of our applications. The quantities
(k) Q . i
are determined by numerical differentiation (i.e. by differentiating an interpolating polynomial). The formulas, given below, were determined by the following criteria:
Furthermore, we picked discretization formulas that use points at equally spaced intervals. This may be changed in future versions of MICROSCOPE.
Denoting the 'discretization parameter' by h > 0, the above criteria give rise to the following formulas:
(0)
Q = F5
i
(1)
Q = (-F1 + F9)/(2*h)
i
(2) 2
Q = (F1 - 2F5 + F9)/h
i
(3) 3
Q = 4*(-F1 + 2F3 - 2F7 + F9)/h
i
(4) 4
Q = 16*(F1 - 4F3 + 6F5 - 4F7 + F9)/h
i
(5) 5
Q = 243*(-F1 + 4F2 - 5F4 + 5F6 - 4F8 + F9)/(2h )
i
(6) 6
Q = 729*(F1 - 6F2 + 15F4 - -20F5 + 15F6 - 6F8 + F9)/h
i
where
F1 = F(P(i) - hd)
2hd
F2 = F(P(i) - ---)
3
hd
F3 = F(P(i) - --)
2
hd
F4 = F(P(i) - --)
3
F5 = F(P(i))
hd
F6 = F(P(i) + --)
3
hd
F7 = F(P(i) + --)
2
2hd
F8 = F(P(i) + ---)
3
F9 = F(P(i) + hd)
The first four of the above formulas can be found in Abramowitz and Stegun, 1968, p. 914, and the last two can be derived by the techniques described there. Higher order formulas can of course also be derived.
It is interesting to note that
(k) Q iconsidered as a function of x (or a function of a continuously varying parameter i) is a function of the same degree of smoothness as the trial function itself. Just the number of points of reduced smoothness is increased, there being one corresponding to each point in the differentiation formula passing through the point of reduced smoothness. Thus numerical differentiation can be thought of as a smoothing process. (This contrasts starkly with its numerical property of a roughing procedure due to round-off effects.) The figures in section 4.2 illustrate the smoothing properties of numerical differentiation.
For a more detailed discussion of the truncation and round-off properties of these formulas see Appendix I.
In this subsection, the alphanumerical screen display generated by MICROSCOPE is described in detail. The purpose of the description is twofold: First, it is necessary for an understanding of the examples sprinkled throughout this manual. Second, however, it also serves as an illustration of some of the capabilities of MICROSCOPE and as a first introduction to the power and versatility of the package.
During a MICROSCOPE session the CRT screen will usually be occupied by the 'alphanumerical display' like the one illustrated in figure 1. (More sophisticated graphical displays can also be obtained, see sections 4.14 and 5.)
222 I : **********************************
. 22 I : * I 22 .
. 222 I : * I 222 .
. 22 I : I 22 .
. 22 I : I 22 .
. 222 I :* I 222 .
.. 22 I : I 22 ..
7654321..876543212287654321.987654321+123456789.123456782212345678..1234567
. 222 I : I 222 .
.. 22 I *: I 22 ..
.. 22 I : I 22 ..
.. 222 I : I 222 ..
... 22I * : I22 ...
.... 22 * : 22 ....
**********************************22222222..........
===========================================================================
CD: deg = 1 dir = ( 0.000000D+00, 1.000000D+00) ch = 6.0000D-04
Point = ( 0.000000D+00, 1.000000D+00) s = 5.0000D-03
Direction = ( 1.000000D+00, 0.000000D+00) h = 3.0000D-02
F0 ( 9.39D-09, 1.27D-02) F2 ( 1.20D-01, 2.22D+00) F3 (-1.20D+01, 1.21D+01)
I/O: 25 6 6 1 2 3 27 28 29 30 NRML on current CALLS = 174
Figure 1 Demonstration of the Alphanumerical Display
In this example, the trial function is f(x,y) = x^2 *y*abs(x*y). We analyze a cross derivative of it in the direction (0,1). The direction of investigation is (1,0), and the point of examination is (0,1). Because of the absolute value term, the cross derivative is twice but not three times tangentially differentiable at the point of examination.
We explain the display starting at the top and referring to lines either by their position or by the first word contained in them.
The display is naturally divided into two parts: The first 15 lines constitute the 'graphical display', and lines 17 through 21 form the 'numerical display'. Line 16 separates the two displays and is present only if there is sufficient space available.
Every MICROSCOPE display shows the graphs of a set of tangential derivatives of degrees between 0 and 6. The function that is being differentiated tangentially is called the 'display function'. It may be the trial function itself (the most frequently occuring case), or, as in this example, a pure cross derivative of a degree between 1 and 6.
Distance along the horizontal axis corresponds to displacement in the direction of investigation from the point of examination (which corresponds to the center of the horizontal axis). The scale on that axis has the same meaning for all graphs shown. The vertical axis corresponds to function values, but the vertical scale is different for each derivative, and is chosen such that the graph fills the entire vertical extent of the graphical display. This normalization is carried out and modified automatically as the investigation proceeds. A period denotes the graph of the display function, and the digit k, say, where k is between 1 and 6 and denotes the k-th tangential derivative. In this example, the display function and its second and third tangential derivatives are being shown.
Sometimes some derivative is more interesting than others, and then it may be useful to accentuate it with asterisks. In the above example, the third derivative has been so emphasized. Usually, the graph of a given derivative overwrites the graph of any lower order derivative. However, the graph of an accentuated derivative overwrites all other graphs.
Both the number of lines and the number of columns in the display can be set according to specific purposes. Throughout this manual we will use 75 columns and 15 rows for the graphical display as this is appropriate for the width of the printed output. A width of 79 might be preferable for a standard CRT terminal displaying 80 columns; 135 columns and a larger number of rows might be used for printed output; and a smaller number is sometimes useful for preliminary investigations when the trial function is expensive to evaluate. An even number of columns is normally undesirable because then the point of examination does not lie in the precise center of the display.
Line 8 in the example contains a horizontal scale counting from the center to the left and right. This will usually be omitted so as not to clutter the display, but it is sometimes useful for the identification of points of interest. The distance from the center can be used to inquire about numerical values of certain derivatives, or to shift the graph left or right as required by the circumstances. The scale can be turned on or off, or replaced with a line of hyphens. The graph of any function overwrites the horizontal scale. The center of the display can optionally be marked with a "+" sign. This is sometimes useful to pinpoint symmetry properties. In the above example the graph of the third derivative passes through the center of the display. The center mark overwrites the graph of any function passing through it.
The central column of colons marks the column corresponding to the point of examination. The two columns of I's outline 'the window', i.e. that set of points at which the numerical differentiation is affected by the behavior of the trial function at the point of examination. (This is not always quite true when cross derivatives are being examined, see the discussion in section 4.11). For every point in the window, the interval about the point spanned by the differentiation stencil contains the point of examination. The differentiation formulas are chosen such that the window is identical for all tangential derivatives. Obviously, for the display function itself the window contains only the point of examination. For its derivatives, the window contains all points that can be written as P + td, where abs(t) does not exceed the value of h/s, (and, as before, P is the point of examination, d is the direction of investigation, h is the discretization parameter, and s is the interval). We define w:=2*h/s to be the 'window width'. Obviously, by picking the window width < 1, it is possible to generate displays in which the only point contained in the window is P itself. However, for several reasons it is usually preferable to use a larger window width. This aspect is discussed in more detail in section 4.1. The default value of the window width is 12, as in the example, but this can be changed as needed.
Without examining the numerical display, it is apparent from the graphical display that some function is being investigated that is twice but not three times differentiable since the (accentuated) third derivative shows a clear step discontinuity. Note how that step is smeared out over the window.
The first line of the numerical display starting with " CD:" indicates that the derivatives are mixed partials, i.e. they are tangential derivatives of a cross derivative (of the trial function). This line is absent if the trial function itself is being examined, and the line carries only information that is pertinent to the cross derivative. In the present example, the degree of the cross derivative ( "deg") is 1, the cross direction is (0,1), and the discretization parameter used for computing the derivative (denoted by "ch" corresponding to "h" for the tangential derivative) is 6D-4. The formulas for computing cross derivatives are equivalent to those for computing tangential derivatives.
The line starting with "Point" describes the point of examination (which is (0,1)) and the interval (which is 5D-3). This display implies that the trial function is bivariate. Appropriate modifications take place for functions of 1 or 3 variables. The next line starting with "direction" describes the direction of investigation (which is (1,0)) and the discretization parameter (which is 3D-2). Notice that the window width is 2h/s = 12 which is consistent with the graphical display
The "Point" and "Direction" descriptions are always present.
The next line, starting with "F0", gives the ranges of the displayed derivatives. Ranges are labeled Fk where k is the order of the tangential derivative. For example, F2 corresponds to the range of the second (tangential) derivative (of a first order cross derivative) which here is (0.12,2.1). If more functions (up to seven) are displayed, then their ranges are also given, adding up to two lines to the numerical display as necessary.
The next and last line, starting with "I/O", contains information about the present configuration of MICROSCOPE.
The 10 integers following "I/O" are the FORTRAN channel numbers as they are employed by MICROSCOPE at the display time. They correspond to the following:
The first six I/O channel numbers are parameters of the master MICROSCOPE routine (see section 3.1) and can each be changed during a session. The last four channels are zero (i.e. unused) by default, and can be set to non-zero (i.e. active) values during a MICROSCOPE session. It is of course unnecessary to have all channel numbers distinct. For example, all three loading devices might be identical, allowing for the possibility of having just one list of points, directions, and cross directions.
Following the device numbers is a flag "NRML on" or "NRML off" (here on), indicating whether the direction of investigation entered into MICROSCOPE has been normalized or not. The "Direction" line always gives the direction as entered by the user. If NRML is off, the direction of investigation, d, is identical to that printed in the display, otherwise it is that printed in the display normalized to have unit Euclidean length. With the normalization on, tangential derivatives are the standard directional derivatives as they are defined in Calculus Textbooks. With the normalization off, they are the more versatile Gateaux derivatives as they are used frequently in multivariate interpolation and approximation. The NRML flag applies similarly to the cross direction.
Following the NRML flag is the word "current". This indicates that the numerical display corresponds to the graphics display. To explain why this might not be so requires a brief excursion into the Organization of MICROSCOPE. When drastic parameter changes are taking place it is often more efficient to accumulate them before generating the new display. Usually in that case the alphanumerical display is scrolled while commands are being given, and vanishes from the screen. However, sometimes it is preferable to view the old display rather than the commands given so far. In that case, the display can be recalled on the screen, and the word "current" is replaced by "GO pndg" (i.e. the GO command needs to be given to incorporate into the graphical display changes of parameters that may already be shown in the numerical display).
The last item in the last line is "CALLS = 174" where 174 is the number of evaluations of the trial function expended up to this point in the MICROSCOPE session. This allows monitoring the computational effort spent. A typical multivariate surface will be expensive to evaluate, and its cost will usually dominate the cost of running MICROSCOPE and thereby determine how long you have to wait for a reponse from the terminal. Therefore, a significant effort has been made to keep the number of necessary function evaluations small. This feature contributes in no small measure to the complexity of the MICROSCOPE source code. The example display has 75 columns and we are examining a first order cross derivative of the trial function. Each evaluation of the cross derivative requires two evaluations of the trial function. Displaying the graph of the cross derivative alone would therefore require a value of CALLS equal to 2*75 = 150. Note that computing and displaying the second and third tangential derivatives of the cross derivative requires an additional effort of only 16%. This is not a trivial accomplishment as may be suggested by the fact that the same display with the slightly different window width of 10 would require 330 evaluations. For a detailed discussion of the numerical effort associated with various window widths see section 4.3.
Notice that in the above example we identified a step discontinuity in a mixed partial derivative of fourth order (a third order tangential derivative of a first order cross derivative) without fuss or doubt. If this was a real investigation we could state now that the underlying scheme in any case is at most three times differentiable.
With MICROSCOPE, one answers questions like those in section 1.1 for specific examples. Thus only negative results (e.g. that a scheme is not differentiable) can be proved strictly. However, positive results can also be obtained. If the examples are chosen carefully to exclude artifacts and make them representative, our program allows for answers to the above questions with a degree of confidence that borders on certainty! Naturally, a computational approach does not replace a rigorous proof. However, in the authors experience it is often possible to confirm (or shatter) fuzzy notions with MICROSCOPE, and to generate further insight into the phenomena under consideration. Encouraging MICROSCOPE results may help in building up the stamina needed for attempting a rigorous proof, and detailed investigations may help in providing ideas for carrying out the analysis.
Often, of course, MICROSCOPE will be employed in hindsight, when the analyst already knows the answers to questions like the above from theoretical reasoning. The examination with MICROSCOPE then serves to corroborate the reasoning, and also, more mundanely but no less importantly, to eliminate bugs in the implementation.
In the authors' experience, MICROSCOPE proved invaluable as a debugging tool, often pinpointing within narrow limits that part of a code that was at fault, but also sometimes uncovering flaws in the theoretical approach, and, on a few exhilarating occasions, leading to the discovery of new and unlooked for results that could then be proved theoretically. Thus MICROSCOPE is useful in the identification as well as in the verification of results.
In the development of MICROSCOPE we were occasionally confronted with the objection that our alphanumerical displays are rather crude. Our primary response to this is that the very purpose of MICROSCOPE consists of translating very subtle properties into displays that make them glaringly apparent. A large part of our interest centers around step functions, and it seems fair to display them as such. Indeed, any curve fitting might disguise the very features in which we are interested. However, we do provide an interface to a sophisticated graphics package. See sections 4.14 and 5 for more information.
In this section, we illustrate the main applications of MICROSCOPE. Consider the simple example of interpolating to a primitive function p by a cubic spline S, say, according to the conditions:
2 3
S(x) = a0 + a1*x + a2*x + a3*x if x < 1
and
2 3
S(x) = b0 + b1*x + b2*x + b3*x otherwise
where the coefficients are defined by
S(0) = p(0), S(1) = p(1), S(2) = p(2), S'(2) = p'(2), S"(0) = 0
and S is twice differentiable at x = 1. Thus we require interpolation to 'position' (i.e. function values) at 0,1,2, and interpolation to the first derivative at 2 by a piecewise cubic function that is twice continuously differentiable. At 0 we impose additionally a "free end" condition. (There is an extensive literature on Splines, see e.g. de Boor, 1978, or Schumaker, 1981)
For our present purposes we assume that p(x) is the exponential, i.e.
x
p(x) = e
where e is the base of the natural logarithm. This function appears to be sufficiently non-polynomial to avoid the introduction of artifacts.
A simple calculation shows that
a0 = 1
2
a1 = (-2*e +12*e-9)/7
a2 = 0
2
a3 = (2*e -5*e+2)/7
and
2
b0 = (5*e -16*e+12)/7
2
b1 = (-17*e +60*e-24)/7
2
b2 = (15*e -48*e+15)/7
2
b3 = (-3*e +11*e-3)/7
In the following examples we will verify that S with the above
coefficients does indeed possess the required properties.
... I : I
... I : I
.... I : I
... I : I
... I : I
.... I : I
.... I : I
.... I + I
..... : I
I ..... I
I :.....I
I : ......
I : I .......
I : I .........
I : I ..........
===========================================================================
Point = 0.000000D+00 s = 5.0000D-03
Direction = 1.000000D+00 h = 3.0000D-02
F0 (-3.33D-02, 6.76D-02)
I/O: 25 6 6 1 2 3 27 28 29 30 NRML on current CALLS = 249
Figure 2 Interpolation to Position at Zero
We notice that the range of the error is from -0.0333 to 0.0676. Since the graph of the error function passes through the vertical axis slightly below the center, this range is consistent with the error being zero, i.e. with interpolation. A more reliable confirmation can be obtained by asking for the value of the displayed function at the center (i.e. 0) of the display. This procedure, which is not shown here, yields a value of -9.4D-11, which is zero up to round-off errors. (If you try to reproduce this example with your own version of MICROSCOPE you may get a slightly different result because of the random nature of round-off errors.) An alternative approach consists of halving the interval, and thereby decreasing the range covered by the graph. In the numerical display the range should shrink correspondingly, maintaining the value of zero in its interior. Halving the interval yields the display:
.... I : I
.... I : I
.... I : I
.... I : I
.... I : I
..... I : I
.... I : I
..... + I
I ..... I
I : .....
I : I.....
I : I ......
I : I ......
I : I .......
I : I .......
===========================================================================
Point = 0.000000D+00 s = 2.5000D-03
Direction = 1.000000D+00 h = 1.5000D-02
F0 (-2.03D-02, 2.88D-02)
I/O: 25 6 6 1 2 3 27 28 29 30 NRML on current CALLS = 287
The range of the error has been reduced to from -0.0203 to 0.0288 and a zero value of the error at the point of examination is still consistent with the graphical display. This procedure can be continued. Proceeding similarly at the point x = 1 yields a similar result, and you might try this for yourself. At the point x = 2, however, we obtain the display:
I : I ..
I : I .
I : I .
I : I .
I : I ..
I : I .
I : I ..
.. I + I .
... I : I ..
.... I : I ..
... I : I ..
.... I : I ...
..... I : I ...
..... I : I ....
......................
===========================================================================
Point = 2.000000D+00 s = 5.0000D-03
Direction = 1.000000D+00 h = 3.0000D-02
F0 ( 2.64D-11, 1.42D-02)
I/O: 25 6 6 1 2 3 27 28 29 30 NRML on current CALLS = 362
Figure 4 Interpolation at x = 2
The graph of the error function touches, but does not cross, the x axis, which is consistent with our interpolating to the first derivative as well as position. The value of the error at the center can be read directly from the numerical display and is -2.4D-11 (you may get a slightly different value because of the random round-off errors). This confirms interpolation to position. Interpolation to the derivative is likely because the graph appears to have a horizontal tangent. However, confidence can be raised by computing the first derivative of the error as well, yielding the display:
I : I 111
I : I 11.
I : I 111 .
I : I 111 .
I : I 111 ..
I : I 111 .
I : I 111 ..
.. I + I 1111 .
... I : I 1111 ..
.... I : 1111 ..
... I :11111I ..
.... I111111 I ...
..... 1111111 : I ...
111111111. I : I ....
1111111111111111 ......................
===========================================================================
Point = 2.000000D+00 s = 5.0000D-03
Direction = 1.000000D+00 h = 3.0000D-02
F0 ( 2.64D-11, 1.42D-02) F1 (-6.11D-02, 1.77D-01)
I/O: 25 6 6 1 2 3 27 28 29 30 NRML on current CALLS = 374
Figure 5 Interpolation to the First Derivative
Inquiring about the value of the first derivative at the center yields the number 5D-4. This is much larger than the value for the error itself, and hence might cause some doubts. However, the derivative is calculated numerically and hence not exactly. The value reported by MICROSCOPE is actually the truncation error. This can be confirmed by halving the discretization parameter yielding the display:
I : I 1111
I : I 111 .
I : I 1111 .
I : I 111 ..
.. I : I 1111 .
.. I : I 1111 .
.. I : I 11111 ..
.. I + I1111 ..
... I : 11111 ..
.. I 11111 I ..
... 11111 : I ..
... 111111 I : I ..
1111111 I : I ....
11111111 .... I : I ....
11111111 ....................
===========================================================================
Point = 2.000000D+00 s = 2.5000D-03
Direction = 1.000000D+00 h = 1.5000D-02
F0 ( 2.64D-11, 3.03D-03) F1 (-4.21D-02, 7.09D-02)
I/O: 25 6 6 1 2 3 27 28 29 30 NRML on current CALLS = 418
Now the value of the derivative is 1.25D-4 which is a reduction from the previous value by a factor 4. This is consistent with the assumption that the reported error actually is the truncation error since it decreases like the square of the discretization parameter h (see the discussion in Appendix I).
The free end condition (S"(2) = 0) can be considered an interpolation condition. However, it cannot be checked by examining the error function, since the value of p"(2) is (in general) unknown. So one has to consider the interpolant itself. We obtain the display:
11 I : I 22222
1 I : I 22222 1
1 I : I 22222 1
11 I : I 22222 11
1 I : I 22222 1
11 I : I 22222 11
11 I : 22222 11
1 I 22+22 I 1
11 22222 : I 11
11 22222 I : I 11
111 22222 I : I 111
222221 I : I 11
22222 111 I : I 111
22222 1111 I : I 1111
22222 1111111111111111111
===========================================================================
Point = 0.000000D+00 s = 2.5000D-03
Direction = 1.000000D+00 h = 1.5000D-02
F0 ( 8.83D-01, 1.12D+00) F1 ( 1.26D+00, 1.27D+00) F2 (-2.53D-01, 2.53D-01)
I/O: 25 6 6 1 2 3 27 28 29 30 NRML on current CALLS = 505
That the second derivative is indeed zero can now be verified using the same techniques as above (i.e. inquiring directly or using a decreased discretization parameter).
In the present example, the only critical set is the one containing the point x = 1 (assuming we are confident that S has been programmed to be polynomial elsewhere). When investigating smoothness it is usually preferable to examine an error (if a primitive function is available) rather than the interpolant itself. In the present example, the error function at the single critical point yields the display:
I : 3333333333333333333333333333333333
I : 3 I .2222
I : 3 I 12222
I : I 112222
I : I .112222
I :3 I .1112222
I : 11122222
I +1112222
I111112222 I
.1111222223: I
1111222222I : I
112222222 I : I
1222222 I 3 : I
1222222 I 3 : I
3333333333333333333333333333333333 : I
===========================================================================
Point = 1.000000D+00 s = 2.5000D-03
Direction = 1.000000D+00 h = 1.5000D-02
F0 ( 2.49D+00, 2.97D+00) F1 ( 2.39D+00, 2.90D+00) F2 ( 2.48D+00, 3.11D+00)
F3 ( 2.73D+00, 4.06D+00)
I/O: 25 6 6 1 2 3 27 28 29 30 NRML on current CALLS = 592
This display can be interpreted similarly as figure 1, and clearly shows a function that is twice but not three times differentiable. This is what we expect from our construction.
The 'precision class' of an approximation scheme is the set of functions that are reproduced exactly by the scheme. Of particular interest are the polynomials in the precision class. The largest number q, say, such that all polynomials of degree up to q are in the precision class is the degree of polynomial precision.
To decide if any particular function is in the precision class one approximates that function and investigates if the error is zero within round-off effects. This applies even to non-polynomial functions.
To answer the more narrow question of polynomial precision it is possible to proceed more systematically and to examine polynomials of increasing degree. The polynomials should be as general as possible, i.e. they should be unlikely to introduce misleading artifacts. One possibility is to choose the coefficients randomly. If really in doubt, one might examine each polynomial in a suitable basis of the appropriate polynomial space. The points of examination should be chosen such that each appropriate part of the domain is covered. For instance, in the present example (which, incidentally, has linear precision) one should choose one point > 1, and one < 1, covering both cases in the definition of S.
The following figure illustrates the typical appearance of a zero error.
. . I .: I .
. . I : I . .
. . I : I ..
. : . I . . .
.I . : I .
I . : I . . .
. . . I . . I.
. . I : I . .
. I : . I . . .
. I :. . I .
. I : I . . .
I : .
. . . . . . . . . I : .I . .
. . I. : I . . . ..
. . . . . I : I . .
===========================================================================
Point = 2.342300D-01 s = 2.5000D-03
Direction = 1.000000D+00 h = 1.5000D-02
F0 (-9.99D-11, 9.78D-11)
I/O: 25 6 6 1 2 3 27 28 29 30 NRML on current CALLS = 667
Notice how the values of the function are spread randomly over the entire display, accurately pinpointing the round-off nature of the error. Sometimes, however, one has a situation in which the coefficients of the interpolant are evaluated inaccurately, leading to an error that has a scattered distribution with a bias like that in the following display:
I : I . . .
I : I . .
I : I . .. .
I : I . . .
I : .I .. . .. . .
I ..: I . . . .
. .I : .. . . . .
. I . :. . I .
.. .. : . I . .
.. . . I . : . I
. .. . . . I . I
. . .. .. I : I
. . .. I : I
. . . I : I
. . . I : I
===========================================================================
Point = 2.342300D-01 s = 2.5000D-03
Direction = 1.000000D+00 h = 1.5000D-02
F0 ( 2.16D-10, 7.31D-10)
I/O: 25 6 6 1 2 3 27 28 29 30 NRML on current CALLS = 742
Figure 10 Round-off With a Bias
Verifying that a certain function is a polynomial of a certain degree is often useful when examining cross derivatives. For example, in the bivariate Clough-Tocher scheme (Alfeld, 1984) which is piecewise cubic on a triangle, differentiability between triangles is forced by the requirement that the normal derivatives across edges be linear (rather than quadratic) along the edges. In this illustration, let us verify that our interpolant is indeed cubic. This can be accomplished by showing that the third derivative is constant or that the fourth derivative is zero (within round-off). One obtains for example the following display:
3 I : I 22222
3 3 I : I 222221
3 I : I 3 3 32222211
3I : I 22222111
3 3 3 I 3 : 3 3I3 222221111 3 3
3 33 3 3 : 3 32322 1111 3 3
3 3 3 3 3 I : 22222. 11113 3 3
33 3 3 I3 23222. 11113 3 3 3 3 3 3 3
22222. 11111 I 3
3 3 3 3 3 22222.I111113 3 3 I 3
22223. 11111 33 :3 I 3 3 3 3 3
3 22223 311111 I : I
22222111111 I : I 3 3 3
222221111 I : I
2222211 3 I : I
===========================================================================
Point = 2.342300D-01 s = 2.5000D-03
Direction = 1.000000D+00 h = 1.5000D-02
F0 ( 1.18D+00, 1.43D+00) F1 ( 1.29D+00, 1.41D+00) F2 ( 3.87D-01, 8.92D-01)
F3 ( 2.73D+00, 2.73D+00)
I/O: 25 6 6 1 2 3 27 28 29 30 NRML on current CALLS = 841
Notice that the third derivative is clearly affected by round-off, without any apparent bias. On the other hand, since the displayed limits of its range are identical the size of the range is less than 1/273th of the derivative value. A closer examination reveals that the third derivative ranges from 2.7306 to 2.7323. A telltale sign of round-off errors is that they can be decreased even further by increasing the discretization parameter. Doubling h yields the display:
I : I 3 322222
3 3 I 3 : I 2222231
33 I : I 22222 111
3 3 33 I : 333I 233223 111
3 3 3 3 : I 3 22322. 31111
3 3 I : I 22222. 111 3
3 3 3 3 I 3 33322222.33 33 1111 3 3 3 3
3 3 3 I 22222. I 1111 3 3
3 33 3 32322. : I1111 3 3 3
22222.I 3 : 1111333 3 3
3 322222. I 11311 I 3
3 3 22222. 11111 : I 3 3
22222. 1111111 I : I
23222. 111111111 I : I
2222211111111 3 I : I 3
===========================================================================
Point = 2.342300D-01 s = 5.0000D-03
Direction = 1.000000D+00 h = 3.0000D-02
F0 ( 1.06D+00, 1.56D+00) F1 ( 1.27D+00, 1.50D+00) F2 ( 1.34D-01, 1.15D+00)
F3 ( 2.73D+00, 2.73D+00)
I/O: 25 6 6 1 2 3 27 28 29 30 NRML on current CALLS = 885
The range is now from 2.731361 to 2.731569 which is a reduction by a factor 8. This is consistent with the third derivative being constant since the round-off errors in the third derivative increase like the third power of h, see Appendix I. Further increases do not alter the basic scattered structure of the third derivative (until the range covered in the display includes the point x = 1, where the third derivative has a step discontinuity). Increasing the value of h should eventually decrease the round-off errors below a level at which the true change in the derivative becomes visible in the display if the third derivative was not constant. The fact that this does not happen indicates that there is no such change.
This section contains a description of the master MICROSCOPE routine MCRSCP and a first introduction to the command language of MICROSCOPE. For a more detailed description consult section 4.
During an interactive session the terminal screen may be either 'active' or 'inactive'. On an active screen, changes in parameters are incorporated immediately into the alphanumerical display which is maintained on the screen. On an inactive screen, commands are accumulated and carried out only when the screen is activated by either the command GO or FORCE. An inactive screen shows the last few commands and the corresponding data whereas an active screen shows at most the current command and its data, in addition to the alphanumerical display. (However, an incomplete alphanumerical display may be called on an inactive screen, see the description of the command RSCREEN in section 4.5 and the discussion of the alphanumerical display in section 1.3.)
A second distinction applies to the way in which screen changes are implemented on an active screen. In one version of MICROSCOPE, which we call the 'screen version', the alphanumerical display is kept stationary and only modified rather than replaced. This version requires you or your system manager to provide four simple screen editing routines that are specific to each computer/terminal combination (see chapter 5.) If this facility is not available, or if you use a hardcopy terminal, then the 'scroll version' of MICROSCOPE may be employed, which scrolls the old screen and replaces it by a new version after each change.
The screen version is much preferable, and in this manual we usually assume it is available.
Another distinction pertains to the utilization of lower case letters. The use of lower case letters in a FORTRAN program is not portable, but now most installations do provide this facility. In this manual we assume lower case letters are available. If so, any lower case letters occuring in the input to the program are treated as if they were upper case. Thus commands can be given in lower case letters. Output from MICROSCOPE consists of a mixture of lower and upper case letters. However, to maintain portability we provide a version of MICROSCOPE that works solely with upper case letters. Section 5 contains more information.
MCRSCP(F, INPUT, OUTPUT, GRAPHC, HELP, RECORD, RESTART,
X LINES, WIDTH, PLOT, NUMRCL, PROMPT)
The first parameter is the name of a DOUBLE PRECISION function (the trial function). It has the form
DOUBLE PRECISION FUNCTION F(X)
DOUBLE PRECISION X
DIMENSION X(1)
.
.
.
F = ...
RETURN
END
and must be declared EXTERNAL in the calling program. MICROSCOPE will set at most the first three components of X before calling F (i.e. it can handle only functions of up to three variables). Notice that F does not have any parameters other than the values of the independent variables. Usually, however, it will also depend on parameters describing e.g. the domain data structure. These are most conveniently passed to F in a common block.
All of the remaining parameters of MCRSCP are input integers. MCRSCP will not change them, and they may therefore be integer values rather than variables. The first six of them (i.e. INPUT, OUTPUT, GRAPHC, HELP, RECORD, RESTART) are FORTRAN device numbers corresponding to the first six I/O numbers listed in the numerical display (see section 1.3). It is your responsibility to set up the necessary correspondences between device numbers and physical devices or disk files.
The last five parameters describe the physical dimensions of the screen display. Their meaning is as follows: (Numbers in parentheses give the value used for the examples in this manual, square brackets contain the allowable ranges for the parameters).
Usually, MCRSCP will be called only once during a session. However, it is possible to call your own routines from MCRSCP. You have to supply a subroutine SUBUSR without any parameters which can be called by MCRSCP A detailed discussion of this procedure is given in section 4.15 on user intervention. For now we may assume that SUBUSR is a dummy routine
SUBROUTINE SUBUSR
RETURN
END
(which is supplied with the MICROSCOPE package and which you can of course write yourself).
When called, MCRSCP first prints an identifying stencil on the OUTPUT device, containing its version number and the last modification date. Next, it checks to see if the screen parameters satisfy the restrictions listed above. If not then MCRSCP exits, rather than terminates, in order to give you a chance to correct your mistake.
On the first call to MCRSCP it computes the machine round-off unit (see section 4.1) and prints its common logarithm to give you an idea of how many decimal digits are carried on the computer. The actual number relevant to your investigation may of course be less depending on the details of the arithmetic involved in evaluating the trial function.
Also on the first call only, MCRSCP then prints the first line of any available news (see section 4.10) to help you decide if you wish to read them.
After setting some defaults MCRSCP then enters the command mode, in which it accepts, interprets and carries out the commands entered on the Input Device.
n >>
where n is the number of the command to be given. It starts at 1 and is incremented by 1 after each command.
A MICROSCOPE command is a string of at least two letters and digits. However, only the first two characters are significant and recognized by the parser. They must be entered in the first two columns immediately following the prompt. On an active screen, this will be in the same line as the prompt, on an inactive screen, or in the scrolling version, on the next line. Blank spaces are significant. (MICROSCOPE simply reads the command using a 2A1 format. The remaining columns can be filled with arbitrary characters. They are not read or processed by MICROSCOPE. This is sometimes useful for including comments in a program. A command line is terminated by pressing the RETURN button.
After a command you will be prompted for the appropriate data (if any are needed). The prompt describes what data are needed and is specific to the command. There are seven types of data requests:
The dimension of the vectors is identical to the number of variables of the trial function. (The default is 2, see section 4.18 for a listing of all defaults). Numbers are processed by an internal parser and are essentially format free. The following distinctions and restrictions apply, however:
After receiving the data, MICROSCOPE processes the command and the data. If the screen is active the command is carried out and the alphanumerical display is updated, otherwise the next command is prompted for.
If a command is not recognized a message is printed to that effect and a new command is prompted for. If data are requested, and a number is not recognized, then new data will be requested. However, if data are syntactically correct, but otherwise meaningless (e.g. you request the graph of the 117th derivative), then the command will be aborted and a new command will be prompted for.
There are currently 63 commands which divide into 11 major groups corresponding to sections 4.5 through 4.15. Rather than describing the groups in this section, we illustrate some of them in the course of a search for a point of reduced smoothness (the 'active point'), giving the names of the commands used, and further illustrating the capabilities of MICROSCOPE.
The trial function underlying this example is
2 x x
f(x) = (eta*x )/2 + e if x > 0 and f(x) = e otherwise
(where eta = -0.005). Clearly f is once but not twice differentiable at the origin. The exponential term serves to conceal the discontinuity in the second derivative. (Similar functions that are 0 through 5 times differentiable can also be examined for arbitrary values of eta.) The parameters are changed by giving the user intervention command USER (assuming the test package has been loaded, see appendix III and section 5).
Let us pretend we do not know the the degree of smoothness but we wish to examine it. We first inform MICROSCOPE that the trial function has one independent variable, using the DMNSN command. Then we enter a search interval (here [-1,3], say) using the IINTVL command. At this stage, the screen is still inactive. This is changed, and the first display is obtained, by giving the GO command. The display will remain active throughout the remainder of this section. We obtain:
I : I ..
I : I .
I : I ..
I : I .
I : I ..
I : I ..
I : I ..
I : I ..
I : I ...
I : I ...
I : I ....
I : I .....
I : .......
............ I
........................... I : I
===========================================================================
Point = 1.000000D+00 s = 5.4054D-02
Direction = 4.000000D+00 h = 3.2432D-01
F0 ( 3.68D-01, 2.01D+01)
I/O: 25 6 6 1 2 3 27 28 29 30 NRML on current CALLS = 960
This display shows no lack of smoothness at all. So, not knowing about the nature of the suspected discontinuity at all we turn on the horizontal scale (using DSCALE) to locate any phenomena of interest and activate (using DGRAPH) the graph of the highest possible (i.e. sixth) derivative. These changes are incorporated one by one (not shown here) into an active screen. The final result is this:
6 I : I ..
I : I .
6 I : I ..
I : I .
6 I : I ..
I : I .666666666
6 6 I : I 666666666666666666
66666666666666321.987654666666666666666666666666123456789.123..6789.1234567
6 6 I : I ...
I : I ...
6 I : I ....
I : I .....
6 I : .......
............ I
....................6...... I : I
===========================================================================
Point = 1.000000D+00 s = 5.4054D-02
Direction = 4.000000D+00 h = 3.2432D-01
F0 ( 3.68D-01, 2.01D+01) F6 (-6.07D+01, 6.27D+01)
I/O: 25 6 6 1 2 3 27 28 29 30 NRML on current CALLS = 972
Now it is clear that there is a phenomenon of some kind occuring at the point -17 measured on the horizontal scale. To examine this further, we shift the graph by 17 (using SHIFT). This centers the interesting parts of the graph, yielding:
I 6 : I ..
I : I .
I 6 : I ..
I : I .
I : 6 I ..
I : I ..
I 6 : 6 I ..6666666666
6666666666666666666666666666666654321.12366666666666666666666666689.1234567
6 :6 I ...
I : I ...
I6 : I ....
I : I .....
I 6: .......
............ I
........................... I 6 I
===========================================================================
Point = 8.108108D-02 s = 5.4054D-02
Direction = 4.000000D+00 h = 3.2432D-01
F0 ( 1.47D-01, 8.00D+00) F6 (-6.07D+01, 6.27D+01)
I/O: 25 6 6 1 2 3 27 28 29 30 NRML on current CALLS = 989
Already at this stage it is possible to tell the nature of the reduced smoothness. The sixth derivative clearly shows four distinct local extrema. This qualifies it as the third (numerical) derivative of a delta function. The second derivative of the trial function must therefore be a step function. The trial function is hence once but not twice differentiable. Derivatives of step functions are discussed in detail in section 4.2 and in appendix I.
To confirm our diagnosis directly we examine the second derivative. Turning on its graph and erasing the distracting graph of the trial function (using EGRAPH), however, yields the following display which shows no lack of smoothness at all in the second derivative:
I 6 : I 22
I : I 2
I 6 : I 22
I : I 2
I : 6 I 22
I : I 22
I 6 : 6 I 226666666666
6666666666666666666666666666666654321.12366666666666666666666666689.1234567
6 :6 I 222
I : I 222
I6 : I 2222
I : I 22222
I 6: 2222222
222222222222 I
222222222222222222222222222 I 6 I
===========================================================================
Point = 8.108108D-02 s = 5.4054D-02
Direction = 4.000000D+00 h = 3.2432D-01
F2 ( 1.48D-01, 8.08D+00) F6 (-6.07D+01, 6.27D+01)
I/O: 25 6 6 1 2 3 27 28 29 30 NRML on current CALLS = 989
The second derivative appears smooth because the discretization is too crude to show its discontinuity. So we keep decreasing the value of the discretization parameter (and simultaneously the interval, keeping the window width fixed). Since the active point is only approximately in the center of the display, this will cause it to drift to one side which has to be counteracted by more shift commands. Because of the organization of MICROSCOPE it is particularly efficient and convenient to divide (or multiply) with powers of two. Halving h four times (using HALVE) followed by a shift of 24 and four more halvings yields, after erasing the scale and the graph of the sixth derivative:
I : I 2222
I : I 2222
I : I 222
I : I 222
I : I 222
222222 : I 222
22 I 2: I 2222
2 222 I 2 I 222
22 2 I :22 I 222
2222 I : 22222
2222 I : I
22 I : I
2222 I : I
222 I : I
2222 I : I
===========================================================================
Point = -1.040834D-17 s = 2.1115D-04
Direction = 4.000000D+00 h = 1.2669D-03
F2 ( 9.92D-01, 1.00D+00)
I/O: 25 6 6 1 2 3 27 28 29 30 NRML on current CALLS = 1353
Now it is clear that there is a jump discontinuity in the second derivative of the trial function. The location of the active point can be seen to be 0 +/- about s = 0.0002. Even the size of the jump can be estimated. Inquiring (using TYPE) for the values of the derivative at the left and right end points of the window (where in this crude display it assumes its extrema) yields values of 9.98846D-1 and 9.96522D-1 respectively, corresponding to a jump of -2.32D-3. The real jump is of course eta = -5.0D-3. The discrepancy is due to a truncation error, not just in approximating the derivative, but also in picking the points of evaluation. The estimate can be improved by extrapolating (linearly) the left and right parts of the second derivative to the point of examination, using an estimate of the third derivative which can be obtained from one of the earlier displays. With the third derivative set equal to 1, this procedure yields an estimate of -2.32D-3 - 12s = -4.86D-3 which is a very good result considering the difficulties of numerical differentiation.
In the present setting, it is not possible to decrease the discretization parameter further without succumbing to round-off errors in the second derivative (you may wish to try this yourself). However, the above example illustrates how to alleviate the limitations of round-off errors. We already recognized the true nature of the reduced smoothness when examining the sixth rather than the second derivative. With the accuracy used for the examples in this manual (10 digits), steps in the second derivative 300 times as small as the one above can be detected. The following display shows the graph of the sixth derivative for the same function as above with eta = -1.7D-5 and h = 0.24. These values have been taken from table 4 in appendix I (tau = 1.0D-10).
I : I 66
I : I 66
I : I 66
I : I 66
I : I 66
I : I 666
I : I 66
I : I 666
I : I 6666
I 6: I 6666
I 6 : 6 I 6666
I : 6 666
666 6 I
666666666666 I 66 : 6 I
666666666666666666 I :6 I
===========================================================================
Point = -1.040834D-17 s = 4.0000D-02
Direction = 4.000000D+00 h = 2.4000D-01
F6 ( 2.23D-01, 4.40D+00)
I/O: 25 6 6 1 2 3 27 28 29 30 NRML on current CALLS = 1440
Figure 18 At the limits of resolution
The sixth derivative shows 4 extrema, identifying the step in the second derivative.
The fact that all arithmetic operations are carried out with finite accuracy is fundamental and all pervasive in working with MICROSCOPE. Numerical Differentiation is particularly sensitive to rounding since it involves the subtraction of nearly equal quantities, leading to a cancellation of significant digits. Consequentially, the error in numerical differentiation does not decrease indefinitely with the discretization parameter, but rather starts to behave erratically beyond a certain minimum value and tends to increase as the discretization parameter is decreased further. However, it is its very randomness that makes round-off errors readily apparent in a MICROSCOPE display so that you will not be misled into accepting meaningless results at their face value.
The precise effect of rounded arithmetic depends on several factors:
Of the above three sources of roundoff errors, you will usually have a firm knowledge only of the first. Knowledge about the other two may be gained by using MICROSCOPE itself. This is usually more effective than attempting a theoretical analysis.
The precise nature of roundoff errors is discussed in Appendix I. Detailed formulas and tables are given. However, the usual approach will be to experiment with various values of the discretization parameter, choosing it as small as possible without generating the telltale random pattern that indicates round-off errors.
One perhaps surprising feature of numerical differentiation is that doubtful discontinuities in some derivatives can often be displayed very clearly by plotting higher order derivatives, even though a large discretization parameter is necessary to keep them free of roundoff errors. This fact becomes quickly apparent empirically, and is analyzed in Appendix I.
The conceptual basis for the approach is provided by a Dirac Delta function (see Courant and Hilbert, 1962, p.456). A step function can be thought of as rising at the discontinuity with infinite slope from one discrete value to another. If it is piecewise constant its derivative can be thought of as a function that is zero everywhere, except at the discontinuity, where it is infinite such that the integral over the derivative equals the jump in the step function. That derivative is a delta function.
More precisely, a delta function about a point P, say, is defined as the limit of a sequence of non-negative functions that are different from zero only in neighborhoods of P. The diameter of those neighborhoods tends to zero, and the value of the integral over the functions remains constant (at 1 in the classical definition).
Derivatives of a delta function can be visualized by considering a function in the above sequence. Suppose it is bell shaped. Then its derivative will have two local extrema of opposite signs. The second derivative will have three extrema, etc. Taken to the limit, the k-th derivative of a step function has k local extrema of alternating signs, all located at the same point. Numerical derivatives of a step function exhibit the same phenomena, if the discontinuity is smeared out over a window containing sufficiently many points. This is illustrated in the pictures below. The step discontinuities occur in the m-th derivatives, where m = 0,1, ..., 5. The discretization parameter is h = 1, and the size of the jump is 1. More precisely, the m-th function is defined by m
s(x) = 0 if x < 0 and x /m! otherwise,
where the exclamation mark denotes the factorial.
Notice that as the degree of the discontinuous derivative increases, the smoothness of the graphs of higher derivatives increases. This is consistent with the interpretation of numerical differentiation as a smoothing process. (The numerical derivatives have the same degree of smoothness as that of the display function, just the number of the points of reduced smoothness is increased.)
1111111..........5555555555444444444466666666665555555555222222222211111111
I : I
I : I
I 6666666666 : I
I : 444444444444444I
I : I
I : 6666666666I
6666666I : 66666666
55555556666666666 : 55555555
444444444444444 : I
I : I
5555555555 : 66666666665555555555I
I : I
I : I
333333333333333222226666666666555555555544444333333333333333........
Point = 0.000000D+00 s = 3.3333D-02
Direction = 1.000000D+00 h = 1.0000D+00
F0 (-1.00D+00, 9.98D-11) F1 (-5.00D-01, 7.68D-11) F2 (-1.00D+00, 1.00D+00)
F3 (-4.00D+00, 4.00D+00) F4 (-4.80D+01, 4.80D+01) F5 (-2.43D+02, 3.65D+02)
F6 (-7.29D+03, 7.29D+03)
I/O: 25 6 6 1 2 3 27 28 29 30 NRML on current CALLS = 1662
Figure 19 Derivatives of a Step Function (m = 0)
222222222211 665 444 66 33333 2222222222
I 221111 65 65 4 : 4 6 36 33 22 I
I 22 1111 665 655 4 : 4 633 66 33 22 I
I 22 11615 54 : 4 33 6 5 33 I
I 22 6 51111 6 45 : 436 6 5255533 I
66666666 266 5 1116 55 : 33 6 662 5553 66666666
I666666 6 25 4 11115:33 4 22 6 666666I
555555555 6666 5 22 4 6 3511 6 4 22 5 6666 555555555
I 555 5 22 4 633 :551611 4 22 5 I
44444444 355 24 33 : 5 11142 5 44444444
I444 3555 5 4 2233 6 : 65 22 4115 444I
I 444 335 4 3322 6 : 6 552 45 1111 444 I
I 44433 4 33 22 : 22 5 54 11444 I
I 4443 433 226:622 5 5 4 444 1111 I
I 4443 262 55 444 11111111111
===========================================================================
Point = 0.000000D+00 s = 3.3333D-02
Direction = 1.000000D+00 h = 1.0000D+00
F1 (-1.00D+00, 4.28D-11) F2 (-1.00D+00, 2.85D-10) F3 (-2.00D+00, 2.00D+00)
F4 (-8.00D+00, 1.60D+01) F5 (-8.10D+01, 8.10D+01) F6 (-1.46D+03, 9.72D+02)
I/O: 25 6 6 1 2 3 27 28 29 30 NRML on current CALLS = 1797
Figure 20 Derivatives of a C^0 Function (m = 1)
333333333333322222 666665555555 4444444 3333333333333
I 33 22222 6 56 : 544 44 33 I
I 33 2222 5 : 45 4 33 I
I 3 62 5 6 : 4 5 44 3 I
I 33 6 252 : 4 5 644 I
I 3 5 22 6: 5 666 6664 I
I 3 6 222:4 6 3 6666 I
6666666666666 3 6 5 62 5 6 3 666666666666
555555555 4666 6 5 4: 22 5 6 3 555555555
I 555 46666666 3 :6 22 3 555 I
I 55 4 35 4 : 222 63 55 I
I 55 44 5 3 4 : 6 2625 55 I
I 5 4 33 4 : 33 222 5 I
I 5 455 344 : 6 33 6 55222225 I
I 55555 4444444 33333333366666 55555222222222222222222
===========================================================================
Point = 0.000000D+00 s = 3.3333D-02
Direction = 1.000000D+00 h = 1.0000D+00
F2 (-1.00D+00, 1.99D-10) F3 (-1.00D+00, 1.74D-09) F4 (-2.67D+00, 2.67D+00)
F5 (-8.98D+00, 1.35D+01) F6 (-1.46D+02, 1.46D+02)
I/O: 25 6 6 1 2 3 27 28 29 30 NRML on current CALLS = 1932
Figure 21 Derivatives of a C^1 Function (m = 2)
444444444444444443333 666 5555555 44444444444444444
I 444 3333 6 : 6 5 5 444 I
I 4 333 6 : 6 5 554 I
I 44 33 6 : 6 445 I
I 4 33 : 5 4 55 I
I 4 633 : 5 6 4 5 I
I 4 33:5 4 555 I
55555555555555 4 6 53 6 4 55555555555555
I 555 4 5: 33 4 I
I 55 4 6 5 : 33 6 4 I
6666666666666666 5 6 5 : 336 6666666666666666
I 66 5 4 5 : 433 66 I
I 66 55 6 45 : 44 6333 66 I
I 6 5 6 55 4 : 4 6 33336 I
I 666666555 4444444 666666333333333333333333333
===========================================================================
Point = 0.000000D+00 s = 3.3333D-02
Direction = 1.000000D+00 h = 1.0000D+00
F3 (-1.00D+00, 1.60D-09) F4 (-1.33D+00, 1.14D-08) F5 (-3.00D+00, 3.00D+00)
F6 (-1.28D+01, 2.70D+01)
I/O: 25 6 6 1 2 3 27 28 29 30 NRML on current CALLS = 2067
Figure 22 Derivatives of a C^2 Function (m = 3)
555555555555555555544444 : 666666 5555555555555555555
I 55 4444 : 6 6 55 I
I 5 44 : 6 66 5 I
I 5 44 : 65 I
I 5 44 : 6 566 I
I 5 4 :6 5 6 I
I 5 44: 5 666 I
66666666666666666 5 6 5 66666666666666666
I 666 5 :44 5 I
I 66 5 6: 4 5 I
I 6 5 6 : 44 5 I
I 6 5 6 : 45 I
I 66 55 : 55 444 I
I 6 66 : 5 444 I
I 66666 5555555 444444444444444444444444
===========================================================================
Point = 0.000000D+00 s = 3.3333D-02
Direction = 1.000000D+00 h = 1.0000D+00
F4 (-1.00D+00, 9.60D-09) F5 (-1.38D+00, 6.61D-08) F6 (-4.13D+00, 4.13D+00)
I/O: 25 6 6 1 2 3 27 28 29 30 NRML on current CALLS = 2202
Figure 23 Derivatives of a C^3 Function (m = 4)
6666666666666666666665555 : 666666666666666666666
I 66 555 : 66 I
I 66 55 : 66 I
I 6 55 : 6 I
I 6 55 : 6 I
I 6 5 : 6 I
I 55: I
I 6 5 6 I
I 6 :55 6 I
I 6 : 5 6 I
I 6 : 556 I
I 6 : 655 I
I 6 : 6 55 I
I 6 : 6 555 I
I 66666 5555555555555555555555555
===========================================================================
Point = 0.000000D+00 s = 3.3333D-02
Direction = 1.000000D+00 h = 1.0000D+00
F5 (-1.00D+00, 5.09D-08) F6 (-1.65D+00, 2.67D-06)
I/O: 25 6 6 1 2 3 27 28 29 30 NRML on current CALLS = 2337
Figure 24 Derivatives of a C^4 Function (m = 5)
Sometimes it is desirable to use very small window widths or even completely separate differentiation stencils, for example when searching a large domain for a point of reduced smoothness. Then steps in a derivative will appear as such rather than as gradual slopes that might be obscured by other features of the relevant graphs. Also, cusps and poles are best displayed using small windows (see appendix III for some examples). Usually, however, it will be preferable to use larger windows.
There are two main reasons to use overlapping differentiation stencils: first, as illustrated in the preceding subsection, this is necessary to resolve all the extrema of a derivative of a delta function, and, second, it allows the use of one evaluation of a trial function for several values of a derivative, thereby increasing efficiency and accelerating your interaction with MICROSCOPE. Using overlapping stencils is equivalent to having a window width (w = 2h/s) that is greater than 1.
Due to the internal organization of MICROSCOPE, for any given window width, the number of evaluations (of the display function) required is determined solely by the largest degree of the plotted derivatives. The following table lists formulas for computing the number N of evaluations required given the window width (w), the highest degree of a plotted derivative (k) and the number of columns (C) in the graphical display.
max k: 0 1 or 2 3 or 4 5 or 6
w=2h/s
< 1 C 3C 5C 9C
2 C C + 2 2C + 3 4C + 5
4 C C + 4 C + 4 3C + 8
6 C C + 6 2C + 9 2C + 9
8 C C + 8 C + 9 3C + 16
10 C C + 10 2C + 15 4C + 25
12 C C + 12 C + 12 C + 12
14 C C + 14 2C + 21 4C + 35
16 C C + 16 C + 16 3C + 32
18 C C + 18 2C + 27 2C + 27
20 C C + 20 C + 20 3C + 40
22 C C + 22 2C + 33 4C + 55
24 C C + 24 C + 24 C + 24
Table: Efficiency of Window Widths
The special role played by w = 12 is obvious: the numerical effort is independent of the degree of the plotted derivatives, and exceeds the number of columns in the display only by a constant. This is true whenever w is a multiple of 12. The reason for this is, of course, that the display function is sampled at multiples of 1/2 and 1/3 of the discretization parameter. For efficiency, the default value of w is 12. This is also sufficient for most examinations involving delta functions. Sometimes, however, a larger multiple, like 24 or 36 may be desirable.
This subsection contains a quick reference table of all available commands. For each command, a one-line description is given. Often, the information contained in these short descriptions will be all that is needed. Therefore, the descriptions are used in several other places An appropriate selection of them occurs at the beginning of each of the following subsections. The HSUMRY command lists them ordered by topics. The more detailed HELP command gives a documentation of an individual command headed by the one-line description.
Name Description Data Error Screen
ACCENT mark derivative with asterisks TOGGLE k -1<k<7 immediate
C1CROSS read 1st component of cross direction C1 C><0 deactiv.
C2CROSS read 2nd component of cross direction C2 C><0;1 deactiv.
C3CROSS read 3rd component of cross direction C3 C><0;1 deactiv.
CCHANNL change channel number ic,nc 4 ---
CDIRCTN read direction of cross direction C C><0 deactiv.
CHVALUE read ch for cross differentiation ch 5 deactiv.
CORDER read order of cross derivative k -1<k<7 deactiv.
CWINDOW change window width by factor SIGN m m=0 immediate
D1CHDIR read 1st component of dir. of inv. D1 D><0 deactiv.
D2CHDIR read 2nd component of dir. of inv. D2 D><0;1 deactiv.
D3CHDIR read 3rd component of dir. of inv. D3 D><0;1 deactiv.
DCENTER draw center of graph. display TOGGLE --- --- immediate
DGRAPH draw graph k -1<k<7 immediate
DIVIDE divide h by integer m SIGN m d><0 immediate
DMNSN read dimension v of domain v 0<v<4 deactiv.
DOUBLE double h --- --- immediate
DSCALE draw scale in graphical display TOGGLE --- --- immediate
DXAXIS draw horizontal axis TOGGLE --- --- immediate
EGRAPH erase graph k -1<k<7 immediate
EXIT return to calling program --- --- deactiv.
FLIP replace d by -d TOGGLE --- --- immediate
FORCE recompute everything --- --- activ
GO compute whatever is needed --- --- activ
HALVE halve h --- --- immediate
HELP obtain detailed help on command cmd cmd 2 ---
HSUMRY obtain a summary of all commands --- --- ---
IDIRCTN read the direction of investigation D D><0 deactiv.
IHVALUE read value of h h h>0 deactiv.
IINTVL read endpoints of line of investig. A,B A><B deactiv.
IPOINT read the new point of examination P --- deactiv.
LCROSS load a new cross direction kc 3 deactiv
LDRCTN load a new direction of investigation kd 3 deactiv
LIST print all currently available commands --- --- ---
LOG give channel number for logging nc --- ---
LPOINT load a new point of examination kp 3 deactiv.
MULTPLY multiply h by factor m SIGN m m><0 immediate
NEWS print the news --- --- ---
NORMAL turn normalization of D on/off TOGGLE --- --- deactiv.
OUTPUT write screen immage on recording device --- --- ---
P1CHPNT read first component of P P1 --- deactiv.
P2CHPNT read second component of P P2 1 deactiv.
P3CHPNT read third component of P P3 1 deactiv.
PAUSE pause until RETURN from channel nc nc --- ---
PLOT enter plotting (<PLOT79>) mode --- --- ---
QUIT terminate MICROSCOPE session --- --- ---
RCROSS read channel number for cross list nc --- ---
RDIRCTN read channel number for direction list nc --- ---
RESTART reestablish earlier parameter settings --- --- ---
ROTATE change direction of investigation DD D><0 deactiv.
RPOINT read channel number for point list nc 1 ---
RSCREEN refresh screen --- --- pending
RWIND rewind device nc nc --- ---
SETDF reset defaults --- --- deactiv.
SHIFT shift point of examination is --- immediate
STORE store current parameter settings --- --- ---
TCENTER type value of k-th derivative at center k -1<k<7 ---
TNOTE take a note on recording device --- --- ---
TYPE Type value of derivative at some point k,is -1<k<7 immediate
UNDO undo previous parameter changes --- --- ---
USER enter user mode --- --- ---
WAIT make program wait for GO or FORCE com. --- --- deactiv.
ZOOM change window width by factor SIGN m m><0 immediate
Quick Reference Table of commands
TOGGLE the command reverses itself if called again
SIGN means that reversing the sign of the input parameter
generates the opposite effect (e.g. multiply instead of
divide).
-- Data --
A one endpoint of the line of investigation
B the other endpoint of the line of investigation (P = (A +
B)/2)
C cross direction
C1 first component of cross direction
C2 second component of cross direction
C3 third component of cross direction
ch discretization parameter for cross derivatives
cmd command name
D direction of investigation
D1 first component of direction of investigation
D2 second component of direction of investigation
D3 third component of direction of investigation
DD correction to be added to D
ic type of device (1 = input, 2 = output, 3 = graphic, 4 =
record, 5 = restart)
is integer shift in number of columns on display, +ve shifts
graph to right, point to left, -ve is reverse
k degree of a derivative (to be plotted, erased, etc.)
kc position of the cross direction in a list
kd position of the direction of investigation in a list
kp position of the point of examination in a list
m integer multiplier.
n integer divider
nc channel (device) number
P point of examination
P1 first component of the point of examination
P2 second component of the point of examination
P3 third component of the point of examination
--- no data are required
Restrictions on the data are listed in this column. The entry "---" means that no check is made for errors. The symbol "><" means "not equal". If the error conditions cannot be given explicitly then the following reference numbers have been used:
activ the command activates the screen
deactiv. the command deactivates the screen - i.e. the GO or FORCE
command must be given before the computation proceeds. For
some commands, there is a check if the newly input datum is
identical with the old one, in which case the screen remains
active.
immediate If the screen is active then the effects of the command are
incorporated before the next command is requested
pending If the screen is inactive, this command will put the
graphical display on the screen, even if computations are
pending.
--- The graphical and numerical display, and the screen status
are not changed. Notice, however, that the actual screen
display may be destroyed by error messages or help
information. The next display changing command (or RSCREEN)
will then refresh the screen.
In this and the following subsections, input to MICROSCOPE is via the assigned input channel, unless otherwise stated. This will usually be the terminal, but may also be a file for a batch type session.
ACCENT mark derivative with asterisks TOGGLE k -1<k<7 immediate DCENTER draw center of graph. display TOGGLE --- --- immediate DGRAPH draw graph k -1<k<7 immediate DMNSN read dimension v of domain v 0<v<4 deactiv. DSCALE draw scale in graphical display TOGGLE --- --- immediate DXAXIS draw horizontal axis TOGGLE --- --- immediate EGRAPH erase graph k -1<k<7 immediate FLIP replace d by -d TOGGLE --- --- immediate RSCREEN refresh screen --- --- pending TCENTER type value of k-th derivative at center k -1<k<7 --- TYPE type value of derivative at some point k,is -1<k<7 immediate
We now turn to the descriptions of individual commands. Screen display commands control what appears on the terminal screen. (For the effects of these commands on the < PLOT79> drawing see section 4.14.) There are several subgroups of commands:
By default, only the display function is displayed. The graph of the display function itself is marked by periods, the graph of its k-th tangential derivative (k = 1,2, ..., 6) by the digit k. Using DGRAPH or EGRAPH respectively, you can turn on or off the graphs of the zeroth through sixth derivatives. It is also possible to display no functions at all. In any case, the display function is evaluated only at the points where it is needed.
If you wish to emphasize the graph of a particular derivative then you can plot its graph using asterisks rather than a period or the digit k. The command ACCENT turns this option on or off. Note that the numerical display contains no direct information on which derivatives have been accentuated. Also, if you accentuate several derivatives at once their graphs cannot be distinguished from each other. The graph of an accented derivative overwrites anything else in the display, except for the center mark. The ACCENT command also has the effect of turning the graph of a derivative on, i.e. it does not need to be preceded by a DGRAPH command. By default no derivative is accented.
A horizontal axis can be drawn through the center of the graphical display, using the command DXAXIS. If desired, a scale can be obtained by using DSCALE. Both commands are on/off switches. The point at which the axis or scale intersects the central column of the graphical display does not in general correspond to the value zero. (That point is different for each of the tangential derivatives plotted, and may not even be contained in the graphical display.) The purpose of the axis is to give help in examining symmetry properties, and the scale is useful in locating points of interest for more detailed examination. Both the scale and the axis are off by default.
The on/off command DCENTER may be used to mark the center of the graphical display with a plus sign. This is sometimes useful for confirming symmetry properties. The center mark is off by default.
The FLIP command has the effect of replacing the direction of investigation d by -d. This turns the graph of any derivative about the vertical axis, and it also turns the graph of any odd degree derivative about the horizontal axis, thereby possibly disentangling graphs that have been overwritten by others.
Consider for example the exponential function and its first derivative. Ordinarily the two graphs would be identical, yielding a display in which the display function is completely covered by the graph of its first derivative:
I : I 111
I : I 111
I : I 1111
I : I 111
I : I 1111
I : I 1111
I : I 1111
I : I 1111
I : 11111
I 11111 I
111111 : I
111111 I : I
1111111 I : I
11111111 I : I
111111111 I : I
===========================================================================
Point = 0.000000D+00 s = 1.6667D-02
Direction = 1.000000D+00 h = 1.0000D-01
F0 ( 5.40D-01, 1.85D+00) F1 ( 5.41D-01, 1.86D+00)
I/O: 25 6 6 1 2 3 27 28 29 30 NRML on current CALLS = 2424
After using the FLIP command one obtains:
... I : I 1111111111
... I : I 1111111
.... I : I 1111111
... I : I1111111
.... I : 11111
.... I 11111 I
.... 11111 : I
11111 I : I
1111 ..... : I
1111 I ..... I
111 I : ......
1111 I : I ......
111 I : I .......
111 I : I ........
111 I : I .........
===========================================================================
Point = 0.000000D+00 s = 1.6667D-02
Direction = -1.000000D+00 h = 1.0000D-01
F0 ( 5.40D-01, 1.85D+00) F1 (-1.86D+00,-5.41D-01)
I/O: 25 6 6 1 2 3 27 28 29 30 NRML on current CALLS = 2424
Figure 26 The Disentangled Graphs
where the two graphs are now clearly visible.
The alphanumerical display does not provide precise numerical values. Values of any computed derivative to 21 digits may be obtained by using the TCENTER or the TYPE command. TCENTER provides values at the center of the display, TYPE does so for any point in the display. For TYPE, measure the distance in columns form the center of the graphical display, positive to the right and negative to the left. For clarity, if the screen version of MICROSCOPE is employed, the column in which the value is read is briefly marked with exclamation marks. Returning to command mode requires typing a carriage return. This is necessary so that the prompt for the next command does not overwrite the numerical information before it can be read.
Only derivative values that have been previously computed can be read. This includes all derivatives through [k] where k is the highest degree of a derivative that has been previously displayed on the current active screen, and [k] is k+1 if k is odd or k if k is even.
The command RSCREEN has two distinct applications:
On an active screen, it serves to refresh the screen image after it has been destroyed e.g. by a system message or by HELP information.
On an inactive screen, the last few commands instead of the alphanumerical display are shown. In that situation, RSCREEN calls the display onto the screen without carrying out any computations. Thus changed parameters may be shown in the numerical display without having been incorporated into the graphical display. This condition is indicated by the flag "GO pndg" in the numerical display. (If the screen is active then the flag "current" is displayed instead). The second application is useful if the alphanumerical display has to be viewed in order to decide on the next command. When entering MCRSCP the screen is inactive.
The DMNSN command informs MICROSCOPE of the number of independent variables (1,2 or 3). It is listed here because it does affect the numerical display. However, this command is usually given at most once at the beginning of the MICROSCOPE session. The default dimension is 2. DMNSN also resets the point of examination, the direction of investigation, the cross direction, and the order of the cross derivative to their default values.
IINTVL read endpoints of line of investig. A,B A><B deactiv. IPOINT read the new point of examination P --- deactiv. LPOINT load a new point of examination kp 3 deactiv. P1CHPNT read first component of P P1 --- deactiv. P2CHPNT read second component of P P2 1 deactiv. P3CHPNT read third component of P P3 1 deactiv. RPOINT read channel number for point list nc --- --- SHIFT shift point of examination is --- immediate
The simplest way of defining the point of examination is to read it from the input device by using the IPOINT command. Individual components may be defined similarly by using the commands P1CHPNT, P2CHPNT or P3CHPNT.
Sometimes it is useful to generate a list of possible points of examination before the MICROSCOPE session, and then refer to them by their position (i.e. line number) in the list, rather than by typing the individual components. This is facilitated by the LPOINT command, which prompts for the position of the point and then reads it from a device (usually a file) with a device number that has been set by the RPOINT command (the default is zero). All point loading devices must contain one vector per line.
The SHIFT command can be used to shift the graph by a specified number of columns to the right (if the number is positive) or left (it is negative). There is no restriction on the number of columns by which to shift the graphs. In particular, the number may exceed the number of columns in the display.
It is also possible to define the part of the domain that is to be shown in the display by defining its two endpoints. This is accomplished by the IINTVL command, which will prompt for the left and right endpoint (in that order) of the graphical display. These have to be separated by at least one line feed. Notice that IINTVL defines simultaneously the point of examination, the direction of investigation, the discretization parameter, and the display interval.
D1CHDIR read 1st component of dir. of inv. D1 D><0 deactiv. D2CHDIR read 2nd component of dir. of inv. D2 D><0;1 deactiv. D3CHDIR read 3rd component of dir. of inv. D3 D><0;1 deactiv. IDIRCTN read the direction of investigation D D><0 deactiv. LDRCTN load a new direction of investigation kd 3 deactiv NORMAL turn normalization of D on/off TOGGLE --- --- deactiv. RDIRCTN read channel number for direction list nc --- --- ROTATE change direction of investigation DD D><0 deactiv.
The direction of investigation can be read from the input device by using the command IDIRCTN. Individual components can be changed by the commands D1CHDIR, D2CHDIR or D3CHDIR.
Similarly to the point of examination, the direction of investigation may also be defined by giving its position in a list. The command LDRCTN loads the direction from the list (usually a file) whose channel number has been defined by the RDIRCTN command.
A new direction may also be defined by rotating the old one. The command ROTATE will read the difference between the old and the new direction, and update the direction. This is sometimes useful for examining the effects of deviating slightly from the previous direction of investigation.
The numerical display will always contain the direction of investigation as it has been input to MICROSCOPE. Internally the direction of the derivative (i.e. the vector d in the differentiation formulas) is usually the given direction divided by its Euclidean length. This leads to the standard definition of a directional derivative. However, sometimes it is more useful to use the unnormalized direction for the differentiation (i.e. to have d identical to the displayed numerical value). The normalization can be turned on or off using the command NORMAL, and is on by default. Normalization changes the values of the directional derivatives, and the domain that is covered by the graphical display. An identical graphical display (but with different numerical values) can be obtained by dividing the display interval instead of the specified direction of investigation by the same factor.
CWINDOW change window width by factor SIGN m m=0 immediate DIVIDE divide h by integer m SIGN m d><0 immediate DOUBLE double h --- --- immediate HALVE halve h --- --- immediate IHVALUE read value of h h h>0 deactiv. MULTPLY multiply h by factor m SIGN m m><0 immediate ZOOM enlarge central portion of displ. SIGN m m><0 immediate
This group of commands controls the three parameters h (discretization parameter), s (the display interval), and w (the window width) which are related by w = 2h/s (w must lie in the interval [1/16,576] and can assume only certain rational values). Each of the commands changes two of the three parameters, leaving the third constant, according to the following table:
Parameter left constant
w s h
DIVIDE CWINDOW ZOOM
DOUBLE
HALVE
IHVALUE
MULTPLY
Table of Discretization Controling Commands
Most frequently, the window width is left unchanged, often throughout the entire MICROSCOPE session. This accounts for the abundance of commands covering that case.
The command IHVALUE serves to define an entirely new value of the discretization parameter (and the display interval).
Often you will want to modify the current value of h, e.g. to get rid of the influence of round-off errors, or perhaps to increase the accuracy of numerical values. Due to the organization of MICROSCOPE it is particularly efficient (i.e. requires few evaluations of the trial function) if h is multiplied or divided by powers of 2. The factor 2 is used so frequently that two specific commands (HALVE and DOUBLE) have been provided for this case, which require no arguments.
Multiplications or divisions with other integer factors can be accomplished by using the commands MULTPLY and DIVIDE respectively. It is possible to reverse the meaning of the commands by entering negative factors (e.g. MULTPLY followed by -3 is equivalent to DIVIDE followed by +3). This is due more to the internal workings of MICROSCOPE than to any attempt to provide a truely meaningful convenience.
The command CWINDOW changes the current window width by the specified factor. A positive factor corresponds to an increase of the window width, a negative one to a decrease. The change is obtained by altering the value of the discretization parameter while leaving the display interval fixed.
The command ZOOM has a similar effect. However, it alters the display interval and leaves the discretization parameter unchanged. This is useful when operating close to the round-off level with a value of h that must not be decreased. On the other hand, if round-off errors are irrelevant then the CWINDOW command should be used because it is more efficent.
The primary purpose of the commands described here (other than IHVALUE is to alter the display while the screen is active. If several commands are given while the screen is inactive only the last one will actually be implemented. If several changes are to be affected before the trial function is evaluated, then this can be accomplished by erasing the graphs of all derivatives, and then activating the screen (using GO or FORCE), and making the desired changes.
For the following illustration it is assumed that the screen is activ derivatives through second order have just been displayed, 79 columns are being carried in the graphical display, and the current window width is 12. The table below gives the additional number of evaluations of the display function needed to accomodate changes by factors 2,3,...,8. An entry "--" indicates that the change is impossible (e.g. repeated applications of the HALVE command can only accomplish changes by powers of 2). Entries in parentheses indicate that the change has led to a non-integer half window width and therefore is unlikely to be employed.
2 3 4 5 6 7 8
DIVIDE 46 91 68 91 76 91 80
MULTPLY 46 60 68 91 76 91 80
DOUBLE 46 -- 92 -- -- -- 138
HALVE 46 -- 92 -- -- -- 138
IHVALUE 91 91 91 91 91 91 19
CWINDOW 12 24 36 48 60 163 84
(increasing)
CWINDOW 0 0 (82) (237) 0 (237) (237)
(decreasing)
ZOOM 52 76 96 139 151 163 154
(increasing)
ZOOM 40 52 (116) (182) 66 (198) (204)
(decreasing)
Table: Relative Efforts
EXIT return to calling program --- --- deactiv. FORCE recompute everything --- --- activ GO compute whatever is needed --- --- activ QUIT terminate MICROSCOPE session --- --- --- RESTART reestablish earlier parameter settings --- --- --- SETDF reset defaults --- --- deactiv. STORE store current parameter settings --- --- --- UNDO undo previous parameter changes --- --- --- WAIT make program wait for GO or FORCE com. --- --- deactiv.
To activate an inactive screen the GO or the FORCE command must be used. GO assumes that previous function evaluations are still valid and FORCE reevaluates the trial function at all needed points regardless of prior evaluations. Ordinarily GO will be used for the sake of efficiency, but the FORCE command is necessary if changes have been applied to the trial function (either after leaving MCRSCP or by using the USER command) that invalidate the earlier evaluations.
Usually, the command QUIT is used to terminate a MICROSCOPE session. That command leads to a simple FORTRAN STOP statement. However, it is possible to leave the routine MCRSCP via the EXIT command (leading to a RETURN statement), and then call MCRSCP again, perhaps after changing some screen parameters. After reentering MCRSCP, these will again be checked for consistency, and the screen will be marked inactive. However, the defaults will not be set, so that it is possible to return to MCRSCP and find everything as it has been left with the exception of what you may have reset yourself between MCRSCP calls. After reentering MCRSCP it may be necessary to use FORCE instead of GO to avoid peculiar displays due to altered screen parameters.
The recommended way of changing problem defining parameters (e.g. selecting a new trial function) is to use the command USER (see section 4.15).
Usually the screen is set inactive by some drastic change of key parameters such as the point of examination or the direction of investigation. It can also be set inactive by using the command WAIT. This option is particularly useful in the scrolling version when it may be convenient to effect several changes before scrolling the display just once, rather than having it scrolled after each individual command.
The state of a MICROSCOPE investigation can be stored using the command STORE, and later restored using the command RESTART. STORE writes the values of all relevant parameters onto the Restart Device and then rewinds that device. Thus, if the command STORE is given several times only the status during the last issue of the command can be restored. If it is necessary to store several data settings, then the restart channel number has to be changed by using the CCHANNL command (see section 4.13), or a different file has to be assigned to it by the command USER.
If desired, some of the default settings can be reestablished by using the SETDF command. More precisely, the effects of SETDF are:
Note that any changes regarding a cross derivative or the normalization of the direction of investigation are not undone.
HELP obtain detailed help on command cmd cmd 2 --- HSUMRY obtain a summary of all commands --- --- --- LIST print all currently available commands --- --- --- NEWS print the news --- --- ---
The quickest way of obtaining information on possible commands is the command LIST which will print a list of all currently available commands on the screen. More detailed information on all commands can be obtained through the HSUMRY command which prints the table in section 4.3 reorganized by topics.
Information on a specific command can be obtained by typing HELP, followed by the name of the command you want help on. MICROSCOPE will print not only a short description of the command, but also the current values of those parameters that are affected by the command. Thus if, for example, you wish to find out what the current number of independent variables is you can request HELP on DMNSN.
The command NEWS will print a text that may change from time to time and will also depend on your particular installation. At the time of this writing the "news " are a very brief primer on MICROSCOPE.
Finally, if you type a question mark to MICROSCOPE, it will invite you to enter a command, and will suggest HSUMRY and LIST as possibilities.
LOG give channel number for logging nc --- --- OUTPUT write screen immage on recording device --- --- --- TNOTE take a note on recording device text --- ---
Results obtained with MICROSCOPE can be recorded permanently on the Record Device. The command TNOTE will copy any text you type on the Input Device until it encounters the string 'EC' (End Comments, both E and C must be capital) in columns 1 and 2. It then returns to command mode. The length of text lines may be up to 80 columns or the number of columns specified for the graphical display, whichever is larger. Trailing blanks are stripped, to make the record file more compact.
The command OUTPUT causes MICROSCOPE to write a copy of the graphic and the numerical display onto the recording device. This manual is generated by using TNOTE and OUTPUT throughout a lengthy string of MICROSCOPE commands.
There is also a logging device which is active if its channel number is different from zero (it is inactive by default). The command LOG sets the logging device number. If the device is active all commands, data, displays, and notes are recorded on it.
C1CROSS read 1st component of cross direction C1 C><0 deactiv. C2CROSS read 2nd component of cross direction C2 C><0;1 deactiv. C3CROSS read 3rd component of cross direction C3 C><0;1 deactiv. CDIRCTN read direction of cross direction C C><0 deactiv. CHVALUE read ch for cross differentiation ch 5 deactiv. CORDER read order of cross derivative k -1<k<7 deactiv. LCROSS load a new cross direction kc 3 deactiv RCROSS read channel number for cross list nc --- ---
Cross derivatives, like all derivatives in MICROSCOPE, are calculated using numerical differentiation. The formulas are analogous to those given for tangential derivatives. In general, the cross stencils will have no points in common. Therefore, one evaluation of a k-th order cross derivative requires k+1 evaluations of the trial function.
By choosing the cross direction to be parallel to the direction of investigation the range of tangential derivatives can be shifted from [0,6] to [k,6+k]. This situation, however, is not treated separately by MICROSCOPE and does not lead to savings in function evaluations.
The commands controlling the cross derivative (if any) are similar to those controling the direction of investigation. By default, no cross derivative is computed. It is activated by reading a non-zero order of the derivative, using CORDER, or by defining a cross direction using CDIRCTN (which reads the cross direction), C1CROSS, C2CROSS, or C3CROSS (which read individual components), or LCROSS which loads a cross direction from a channel defined by RCROSS (the default channel number is zero). The default order is 1 if the cross derivative is activated by setting a cross direction. If it is activated by setting the order its default direction is the vector with all of its components equal to 1.
The cross discretization parameter ch is either set explicitly using CHVALUE, or by default, in which case it assumes the same value as the tangential parameter h at the time when the cross derivative is activated. For efficiency, later changes of h have no effect on the value of the cross parameter. Note that any change of ch requires a reevaluation of the display function at all needed points.
Cross derivatives are turned off by setting their order to zero (using CORDER). The NORMAL command has the same effect on cross derivatives as it does on tangential derivatives.
When plotting cross derivatives it is important to keep in mind that the significance of the window may be blurred if a stencil used for the cross differentiation, and centered outside of the window, nonetheless crosses the front. This can be avoided either by keeping the cross discretization parameter much smaller than the display interval, or by choosing the cross direction parallel to the front. In the following example, the trial function is
2
f(x,y) = x *abs(x)*y*abs(y)
the point of examination is [0,1], the direction of
investigation is [1,0], and the cross direction is [1,-1]. The
value of h = 3D-3 and ch = 6D-3. Obviously the effects of the
discontinuity in the second tangential derivative are carried
beyond the window.
111 I33333333333I 22222222222222222222222222222
11 3 : I22 ..111
11 I : 3 ...11
11 3I : 22I3 .... 11
11 I : 2 I .... 11
11 I : 2 I ..... 111
11 3 I :2 I.3....... 11
111 .........2...... 11
11 .......3 I 2: I 3 11
....11 I 22 : I 111
.... 11 3 I 2 : I 11
... 11 I2 : I 3 11
... 311 2 : I 113
... 3 221 : 111
33333333333333333333333332222 I11111111111I 33333333333333333333333333
CD: deg = 1 dir = ( 1.000000D+00, -1.000000D+00) ch = 6.0000D-03
Point = ( 0.000000D+00, 1.000000D+00) s = 5.0000D-04
Direction = ( 1.000000D+00, 0.000000D+00) h = 3.0000D-03
F0 (-7.49D-04, 7.28D-04) F1 ( 1.05D-02, 8.00D-02) F2 (-4.40D+00, 4.18D+00)
F3 (-1.34D+01, 1.00D+03)
I/O: 25 6 6 1 2 3 27 28 29 30 NRML on current CALLS = 2598
Figure 27 An Extended Discontinuity
CCHANNL change channel number ic,nc 4 --- PAUSE pause until RETURN from channel nc nc --- --- RWIND rewind device nc nc --- ---
Some rudimentary programming of MICROSCOPE is possible. The basic approach consists of writing a sequence of commands into a file and assigning it to the input device. The various devices can be switched by using the CCHANNL command. It will prompt for two integers, the first defining the type of channel to be switched, the second giving the new channel number. The types are
The loading devices for the point of examination, the direction of investigation, and the cross direction can be switched by using the commands RPOINT, RDIRCTN and RCROSS, respectively. The Help Device cannot be switched as the information contained in it is an integral part of the MICROSCOPE package and is only read during a session. Indeed, it is advisable to protect the help file against an accidental overwrite if possible, e.g. by using a READONLY switch in an OPEN statement.
Any channel can be rewound by the RWIND command. This facility can be used e.g. to regenerate a specific parameter configuration or to overwrite previously generated information that is no longer needed. An infinite loop can be created by putting a rewind command at the end of the input file.
The PAUSE command can be used to view the output of a MICROSCOPE program at leisure during run time. It will ask for a channel number and then suspend computation until a dummy character is received from that channel. The channel will usually be the terminal.
Computation by a program can be suspended by changing the input device to the terminal (entering interactive mode) and later continued by switching the input channel back to the appropriate device.
Note that for all of the device handling commands no precaution has been taken to ensure that only legal channel numbers are given, and that appropriate devices are assigned to the channel numbers. Comments can be included in a program in two ways. Firstly, a 'comment line' replaces a command line and is marked by a "+" character in column 1 or 2. Secondly, any characters contained in the space starting with column 3 of a command line are ignored by MICROSCOPE. This facility does not exist for lines containing data. Indeed, any spurious characters may lead to the rejection of the data by the parser.
PLOT enter plotting (<PLOT79>) mode --- --- ---
The alphanumerical displays usually generated by MICROSCOPE are adequate to obtain the desired answers to your questions, but they may not be suitable e.g. for publications, reports, etc. Depending upon the availability of suitable hard- and software (<PLOT79>, see Beebe, 1979) you may be able to generate more presentable plots. This is accomplished by the command PLOT which generates a drawing of the current alphanumerical display, augmented by some special options.
Details about the link between MICROSCOPE and <PLOT79> are given in section 5 and your system manager may have put out suitable local information in this respect. For the remainder of this subsection, we will assume that the < PLOT79> package is available, that you have linked it with MICROSCOPE and your own program, and that you have some suitable plotting device available. Appendix IV containes copies of black and white pen plots corresponding to the alphanumerical plots scattered throughout this manual.
The command PLOT puts you into the plot submode. The precise kind of plot you will eventually generate depends upon the MICROSCOPE commands given so far, and the specifications you are going to make in plot mode.
We first discuss the PLOT specifications. After entering the PLOT submode you will see a list of current options on the screen. You can change these options one by one. After each change an updated list will be displayed. To enter an option change enter the reference number of the option followed by the appropriate action (as prompted by the program), to activate the actual plotting process enter 0, and to leave the plotting mode and return to command mode without plotting anything enter any negative number. The list of PLOT options is left unchanged between calls to PLOT.
There are 12 PLOT options:
A title printed above the alphanumerical display, consisting of up to 72 characters can be drawn.
A string of text printed below the alphanumerical display, consisting of up to 72 characters can also be drawn.
For both the title and the legend, a default character height is provided which will be decreased automatically if insufficient space is available.
If ON, labels are drawn giving the degree of the derivative of each graph. On a crowded display, these are sometimes awkwardly placed, and marking the graphs (option 4) may be preferable.
If ON, each point at which a derivative has been evaluated is marked with a special symbol, similarly to the graphical display on the termina The symbols for the k-th derivative are
k: 0 1 2 3 4 5 6 Symbol: dot plus asterisk circle x square triangle
Marking the graphs instead of labeling them has the advantage that it is possible to tell from any small part of the graph to which derivative it corresponds, making it unnecessary to trace a graph from the point of interest to its label. It has the disadvantage that it slows down significantly the speed with which the graphs can be drawn. Also, accentuating a graph is only possible with the marking on (see below). The marking symbol is repeated in the numerical display if the marking option is on, enabling easy association of marks and degrees of tangential derivatives.
The graph of the tangential derivatives can each be drawn in a different color, if supported by the existing hardware. In the most extreme case, when all derivatives are to be plotted in different colors, seven pens for their graphs plus an eighth pen for the numerical display, scales, axes, etc. must be supplied.
MICROSCOPE selects a pen by calling a routine SETPEN(I) where I = 1 corresponds to the pen drawing the axes etc., and I > 1 corresponds to the graph of the (I-2)th tangential derivative (I = 2,3,...,8). This range of I has to be translated into the appropriate <PLOT79> range. On some installations the pens are numbered 0 through 7, on others from 1 to 8. The appropriate SETPEN routines are contained in the files PEN07 and PEN18, respectively (see section 5). If only a smaller range of pens is available you may want to write your own pen selection routine, tailoring it perhaps to the specific range of derivatives in which you are interested. Your routine must be of the following form:
SUBROUTINE SETPEN(I)
INTEGER I,J
C I is in the range from 1 to 8
J = some function of I
C J is in the range of pens available
C on the local plotting device
CALL STPEN(J)
RETURN
END
The subroutine STPEN is a MICROSCOPE routine that calls the &\lt;PLOT79&\gt; routine SETPN(J) if the value of J has changed since the last call to STPEN. SETPEN is not referenced if the color option is turned off.
In the numerical display, the ranges of the tangential derivatives are drawn in the same colors as their graphs.
Using color makes for very nice drawings, but for reproduction purposes it may be preferable to use black and white (like the <PLOT79> drawings in Appendix IV).
A large and prominent integer can be drawn below the legend. We have found this useful, for example, for associating a drawing with the reference number of a point in a list (e.g. the point of examination), and associating or defining a figure number (like in this manual).
No page number is drawn if its value is zero.
A label of up to 20 characters of small type can be drawn in the left bottom corner of the plot.
A similar label can be drawn in the right bottom corner of the plot.
The time at which the plot was generated can be drawn in the top left corner of the plot.
Similarly, the date of the plot generation can be drawn in the top right corner.
Both time and date are printed in the same size as the bottom labels. They are sometimes useful for identification and record keeping.
If desired, the numerical display can be turned off. This is sometimes useful to obtain quick drawings.
If desired, the frame that is normally drawn around the entire <PLOT79> drawing can be turned off.
We now turn to the MICROSCOPE options. Essentially everything that appears in the alphanumerical screen display is also contained in the <PLOT79> drawing, although it may assume a slightly different form.
If a horizontal scale is activated (by the DSCALE command) then a standard labeled horizontal scale is drawn on the bottom of the graphical <PLOT79> drawing. The labels on the scale give the distance from the point of examination in units of the direction of investigation. Notice that this may be the normalized version of the direction given in the numerical display, depending on whether the normalization option is on. If it is, then the scale in the <PLOT79> drawing gives the Euclidean distance from the point of examination.
A horizontal line is drawn through the center of the graphical display if the x axis is turned on (using DXAXIS). This is sometimes useful for illustrating symmetry properties.
If this option is turned on (using DCENTER) the center of the graphical display is marked with a cross.
If a graph is to be accentuated its marks (as discussed above) are drawn twice their normal size. However, if color is available, a better emphasis can usually be obtained by using a special color.
USER enter user mode --- --- ---
Upon receiving the command USER, MICROSCOPE calls a subroutine SUBUSR which is often a dummy routine, but which can be written to suit special purposes. SUBUSR has no parameters, but it can communicate with MICROSCOPE routines or your own routines via COMMON blocks. In utilizing this facility it is important to keep in mind that ordinarily MICROSCOPE has no means of determining that the trial function has changed and that previous evaluations of it are no longer valid. Therefore the command FORCE may have to be used rather than GO.
Five examples of possible applications follow:
If derivatives other than pure tangential derivatives of pure cross derivatives are to be examined this can be accomplished by defining the trial function to be the derivative that is to be displayed, or to be further differentiated, and make it dependent upon parameters that are modified by SUBUSR as needed. In this context, it may be convenient to use the routine
DOUBLE PRECISION FUNCTION DERIV(F,NDIM,P,IORDER,DIR,H)
DOUBLE PRECISION F,P(1),DIR(1),H
INTEGER NDIM,IORDER
...
RETURN
END
which computes at the point P the IORDERth derivative in the direction DIR of the double precision function F, using the formulas given in section 1.2 with the discretization parameter H. F must be declared EXTERNAL in the calling program, and is of the same form as the trial function described in section 3.1. NDIM is the number of independent variables F depends on and must not exceed 3. DERIV is used within MICROSCOPE (to compute any cross derivatives) and can therefore be readily called by your own programs.
Obviously, the modified trial function may itself be a tangential derivative. Then the range of derivatives that can be displayed is shifted from the normal 0th through 6th order to a higher one.
If several trial functions are to be investigated in a single session (or with a single linking and loading process) they can all be written into one DOUBLE PRECISION FUNCTION, and distinguished by a parameter that is passed in a common block, and that may be altered by SUBUSR. The examples in this manual are handled in this manner (see appendix III for details).
Often the trial function depends upon data or parameters that can be changed by SUBUSR without having to leave MCRSCP.
File handling is system dependent. You may find it convenient to handle the assignments of files or devices to channel numbers via USER.
In this subsection, E denotes a vector of 1 to 3 components that are all equal to 1. The default settings for MICROSCOPE are as follows:
No scale or x axis is drawn, the center of the display is unmarked, and no derivative is accentuated. PLOT options: see section 4.14
In this section, which is directed mainly to the local system manager, we describe the files supplied with the MICROSCOPE package, and how they can be put together. Usually, it will be the best procedure to select or write files that generate the most powerful and convenient local version of MICROSCOPE, to put them into one file or a library and to make the appropriate information available to your local users.
MICROSCOPE is a portable set of FORTRAN routines. The user must supply a main program and one or two subroutines that are specific to the application. Thus users must have access to at least the relocatable versions of a suitable subset of MICROSCOPE, if not the source codes.
You should have obtained a magnetic tape containing 21 files in the following sequence.
Html-Note: you can download any of these files by clicking on them. The upper case files are still listed for completeness and consistency, but they are not available for downloading.
Html-Note: The file inter.f contains versions of the screen editing subroutines listed below that work on a Silicon Graphics machine and many other Unix machines. These reoutines were written by Chris Alfeld
Numbers in parentheses give the number of lines (of up to 80 characters each) in the files.
Below, the files are divided into several groups. One or two files have to be picked from each of the groups 1 through 6. Within a group, files are listed by decreasing convenience, and increasing portability. If the bottom files are chosen from each group a completely portable version of MICROSCOPE (except for the use of a few non-FORTRAN characters, namely ?,<,>,!,:,;) is obtained.
LC: Contains those subroutines that require the support of lower case letters in order to work properly. This feature is not portable, but widely available, and greatly increases the convenience in using MICROSCOPE.
UC: To be used in case of LC, if lower case letters are not available. UC is supplied as part of the package, but it can also be generated from LC by replacing all lower case letters with their upper case equivalents, and by replacing the subroutine LCUC by a dummy routine. (Upper case versions of all other files can be generated from their lower case counter parts by replacing lower case letters by uppercase ones everywhere. The leading comments describing the file will of course no longer apply.)
If UC is used in case of LC all output will be in upper case letters, and lower case input will not be recognized.
LC and UC contain their respective versions of the master routine MCRSCP, which is the first routines in each of the two files. All other routines then follow alphabetically.
POS: This file contains a
SUBROUTINE POS(IDEV,N)
which positions the device IDEV at the N-th record, such that subsequent reading from IDEV yields the N-th record. This is accomplished by rewinding IDEV and then reading a dummy character from the first N-1 lines (records) in IDEV. If this at all time consuming it may be worthwhile to replace POS by a faster system dependent routine.
SCREEN: To be used if screen editing is locally available, and highly recommended over its alternative, SCROLL, described below. To take advantage of screen editing, the following four small interface routines must be written for each possible computer/terminal combination (in all cases, IDEVCE is the channel number of the relevant terminal):
SUBROUTINE CLSCRN(IDEVCE,I,J)
INTEGER IDEVCE,I,J
C Clear the screen of the device with number IDEVCE below the point with
C the coordinates (I,J) where I is the horizontal coordinate measured
C from the left of the screen starting at 1, and J is the vertical
C coordinate measured from the top of the screen starting at 1.
...
RETURN
END
SUBROUTINE PCURSR(IDEVCE,I,J)
INTEGER IDEVCE,I,J
C Position the cursor at the position (I,J) of the device with number
C IDEVCE, where I,J are as for CLRSCRN
...
RETURN
END
SUBROUTINE PLCHRS(IDEVCE,I,J,N,CHARS)
INTEGER IDEVCE,I,J,N,CHARS(135)
C Write the string of N characters contained in A1 format in the integer
C array CHARS, starting at the location (I,J). I,J, and IDEVCE are as
C in CLRSCRN.
...
RETURN
END
SUBROUTINE BLSCRN(IDEVCE)
C Blank the screen of the device with number IDEVCE
...
RETURN
END
SCROLL: SCREEN can be replaced with SCROLL if the SCREEN package cannot be used. This will have the consequence that the entire screen display is regenerated after each change, rather than just modified.
PLTLC: This file contains the interface between MICROSCOPE and <PLOT79>. The routines use and recognize lower case letters.
PLTUC: The equivalent of PLTLC, but without the lower case letters. Note, however, that <PLOT79> supports Hershey Fonts that include lower case letters regardless of the computing system. Hence the <PLOT79> drawings generated by PLTUC do contain lower case letters. The communication with the terminal via the command PLOT is completely in terms of upper case letters, however.
With both the PLTLC or PLTUC package a pen selection routine also has to be supplied (see section 4.14). The MICROSCOPE package contains three routines in separate files: PEN18 corresponding to 8 pens labeled 1 to 8, PEN07 corresponding to 8 pens labeled 0 to 7, and PENDMY to be used if color is not available (e.g. on dot matrix printer). However, you may wish to leave it up to the user to select or write his own pen selection routines.
PLTDMY: To be used if <PLOT79> is not available. In this case only alphanumerical screen displays like those shown throughout this manual can be obtained. However, the dummy routine simulates the interaction with <PLOT79> up to the point at which the actual drawing takes place. This feature is useful for the development of MICROSCOPE programs, and makes it possible to run programs without < PLOT79> even if they call the package.
There is only the upper case version of PLTDMY.
Ordinarily, users will not be interested in the specific routine SUBUSR employed in the generation of this manual. They will want to write their own. However, the manual version of SUBUSR, and the routines called by it, are contained in the files:
USRLC: contains the lower case version.
USRUC: contains the upper case version, to be used if lower case letters are not available.
USRDMY: contains a dummy SUBROUTINE SUBUSR, which disables the USER command. If the MICROSCOPE routines are held in a library it may be possible to keep this dummy routine in it as a standby if the user does not supply his own version.
Obviously the <PLOT79> package must be linked with MICROSCOPE if it is to be utilized. This is a large graphics package that can be easily interfaced with almost any graphics device. It can be obtained for a nominal charge from Dr. Nelson H.F. Beebe, Department of Physics, University of Utah, Salt Lake City, Utah 84112, Tel.: 801-581-5254.
The file MAIN contains the driver program as it was used in the generation of this manual. It is listed as a data file because the program needs to be modified to suit local circumstances. In particular, the OPEN statement is not portable.
The package includes two help files that contain relevant help information (but no code). Again, there are the lower and upper case versions:
HELPLC
HELPUC
One of these files (preferably HELPLC) must be readable by the user programs, and the users must be informed how to access them. It is possible to design a sequence of MICROSCOPE commands that overwrites the help file, and, if possible, the user should be denied write access to it. You may also want to keep easily accessible backup copies of these files.
The help files each contain the following items:
MANUAL: containes a verbatim machine readable copy of this manual. It is written with FORTRAN carriage control and ready for printing (at 60 lines per page and up to 75 characters per line). Underlining is accomplished by having the underlines in a separate line with a "+" character in the first column, indicating that the line is to overwrite the preceding one. There is no upper case version of this file. If lower case letters are not available a xeroxed copy of the manual rather than a printed version should be distributed.
This installation guide is chapter five of the manual.
MANLC: This file contains a version of the manual that is suitable for on-line reading and editing. Page breaks and underlines have neen eliminated, otherwise the text is identical to that in MANUAL.
MANUC: The upper case equivalent of MANLC.
The following is a quick key to selecting the appropriate subset of the above files. Proceed by answering the questions and following instructions until you reach a line giving the relevant file names or descriptions. The pen selection routine is described in section 4.14 (and can be supplied by the user if necessary), all other ingredients are described in this installation guide.
1. Are lower case letters available on your machine? Yes GO TO 2
No GO TO 3
2. Have you supplied an interface to screen editing? Yes GO TO 4
No GO TO 5
3. Have you supplied an interface to screen editing? Yes GO TO 6
No GO TO 7
4. Is <PLOT79> available? Yes GO TO 8
No GO TO 9
5. Is <PLOT79> available? Yes GO TO 10
No GO TO 11
6. Is <PLOT79> available? Yes GO TO 12
No GO TO 13
7. Is <PLOT79> available? Yes GO TO 14
No GO TO 15
8. LC, SUPPRT, POS or local equivalent, SCREEN, PLTLC, pen selection
routine, possibly USRLC, screen interface, HELPLC, MANUAL, MANLC
9. LC, SUPPRT, POS or local equivalent, SCREEN, PLTDMY, possibly USRLC,
screen interface, HELPLC, MANUAL, MANLC
10. LC, SUPPRT, POS or local equivalent, SCROLL, PLTLC, pen selection
routine, possibly USRLC, HELPLC, MANUAL, MANLC
11. LC, SUPPRT, POS or local equivalent, SCROLL, PLTDMY, possibly USRLC,
HELPLC, MANUAL, MANLC
12. UC, SUPPRT, POS or local equivalent, SCREEN, PLTUC, pen selection
routine, possibly USRUC, screen interface, HELPUC, MANUC
13. UC, SUPPRT, POS or local equivalent, SCREEN, PLTDMY, possibly USRUC,
screen interface, HELPUC, MANUC
14. UC, SUPPRT, POS or local equivalent, SCROLL, PLTUC, pen selection
routine, possibly USRUC, HELPUC, MANUC
15. UC, SUPPRT, POS or local equivalent, SCROLL, PLTDMY, possibly USRUC,
HELPUC, MANUC
In this appendix, we will discuss in some detail the properties of the formulas listed in section 1.2. Without loss of generality we assume the point of examination is the origin. We also restrict our attention to purely tangential derivatives. Thus we may assume there is only one independent variable. Of particular interest are step discontinuities in certain derivatives (for other types of reduced smoothness see the examples at the end of Appendix III). Thus we consider the trial function
F(x) = PHI(x) + S(m,eta,x)
where PHI(x) is a smooth function (i.e. all its relevant
derivatives are continuous) and S is determined by the
requirements that its m-th derivative be piecewise constant and
have a jump discontinuity of magnitude eta at x = 0.
Explicitly, S is defined by
m
S(m,eta,x) = 0 if x < 0 and eta*x /m! otherwise
where the exclamation mark denotes the factorial.
To simplify matters we will assume S is evaluated exactly, and all round-off errors occur in PHI. The properties of the numerical differentiation formulas depend very critically on the relative sizes of the derivatives of PHI, about which little can be assumed in general. Our approach will be to derive general formulas incorporating the derivatives, but to illustrate the relevant phenomena by using
x
PHI(x) = e
since the derivatives of the exponential are all identical
(assuming the normalization option is on).
We will denote the degree of a derivative by parenthesized superscripts.
A simple calculation shows that, for k >= m, the approximation
(k) (k)
Q of f at x = 0 is given by the fundamental equation
(k) (k)
Q = PHI true value
2 (k+2)
+ h *gamma *PHI truncation error
k
(1)
-k
+ tau*h *(1+abs(PHI))*alpha round-off term
k
m-k
+ eta*delta *h delta term
m,k
If k < m the delta term is absent. We discuss the individual terms in turn:
(k) PHI is evaluated at zero.
This term can be obtained by expanding the error into a Taylor series about h = 0, as is done in any textbook on numerical analysis. The function
(k+2) PHIis to be thought of as having been evaluated at a suitable (unknown) point in the domain of the display. The gamma's are constants given in table 1.
We assume rounding is taking place as described in section 0, i.e. if x is the quantity to be rounded to D digits, say, then x is replaced by
-D -D
z := (1+eps *10 )*x + eps *10 (2)
1 2
where eps and eps
1 2
are random numbers between -1 and +1. The quantity tau in (1)
is a suitable average of the
eps and eps terms in (2). We will 1 2set it to 10^(-D). The function PHI is again evaluated at a suitable point in the display, and the alpha's are the sum of the absolute values of the coefficients in the k-th differentiation formula. Their values are listed in table 1.
We use the expression "delta term" because it arises from differentiating the step function s(m,eta,x), and the first derivative of a step function is a Dirac delta function. The
delta
m,k
are constants that can be read off the MICROSCOPE plots in
section 4.2. They are also given in table 1.
k alpha gamma delta
k k m,k
m = 0 1 2 3 4 5 6
1 2 1/6 1/2 1 - - - - -
2 4 1/12 1 1 1 - - - -
3 24 1/16 4 2 1 1 - - -
4 256 1/24 48 16 3 1 1 - -
5 2430 1/27 365 81 14 3 1 1 -
6 46656 1/36 7290 1458 146 27 4 2 1
Table 1: Constants for Differentiation Formulas
We now address three questions:
We are concerned with the question of when round-off effects begin to obscure the graph of a derivative. More precisely, we ask when the round-off term becomes larger than a certain fraction of the rise in the graph of the derivative
(k) PHIover an interval of length s. An equation for the critical step size, which constitutes a limit on h below which no clear graphs can be obtained, is generated by ignoring the delta term and the truncation error, and approximating the rise of
(k) PHIby the mean value theorem. This leads to:
-k (k+1)
psi *h *(1+abs(PHI))*tau*alpha = 2*h*abs(PHI )/w
1 k
where w = 2h/s is the window width,
(k+1) PHIagain is evaluated at a suitable point, and
psi 1is a safety factor. The larger
psi , 1the clearer will be the graph of
(k) PHI .
Solving for h yields:
(k+1) (1/(k+1))
h = [psi (1+abs(PHI))*tau*alpha *w /(2*abs(PHI ))]
1 k
Table 2 contains threshold values for various values of D (the number of digits carried) and
psi = 5, 1for the case that PHI is the exponential. All derivatives have been evaluated at x = 0.
Digits (D) Degree of Derivative (k)
1 2 3 4 5 6
1 3.46D+00 2.88D+00 3.46D+00 4.34D+00 4.94D+00 6.00D+00
2 1.10D+00 1.34D+00 1.95D+00 2.74D+00 3.37D+00 4.32D+00
3 3.46D-01 6.21D-01 1.10D+00 1.73D+00 2.29D+00 3.11D+00
4 1.10D-01 2.88D-01 6.16D-01 1.09D+00 1.56D+00 2.24D+00
5 3.46D-02 1.34D-01 3.46D-01 6.88D-01 1.06D+00 1.61D+00
6 1.10D-02 6.21D-02 1.95D-01 4.34D-01 7.25D-01 1.16D+00
7 3.46D-03 2.88D-02 1.10D-01 2.74D-01 4.94D-01 8.34D-01
8 1.10D-03 1.34D-02 6.16D-02 1.73D-01 3.37D-01 6.00D-01
9 3.46D-04 6.21D-03 3.46D-02 1.09D-01 2.29D-01 4.32D-01
10 1.10D-04 2.88D-03 1.95D-02 6.88D-02 1.56D-01 3.11D-01
11 3.46D-05 1.34D-03 1.10D-02 4.34D-02 1.06D-01 2.24D-01
12 1.10D-05 6.21D-04 6.16D-03 2.74D-02 7.25D-02 1.61D-01
13 3.46D-06 2.88D-04 3.46D-03 1.73D-02 4.94D-02 1.16D-01
14 1.10D-06 1.34D-04 1.95D-03 1.09D-02 3.37D-02 8.34D-02
15 3.46D-07 6.21D-05 1.10D-03 6.88D-03 2.29D-02 6.00D-02
16 1.10D-07 2.88D-05 6.16D-04 4.34D-03 1.56D-02 4.32D-02
17 3.46D-08 1.34D-05 3.46D-04 2.74D-03 1.06D-02 3.11D-02
18 1.10D-08 6.21D-06 1.95D-04 1.73D-03 7.25D-03 2.24D-02
19 3.46D-09 2.88D-06 1.10D-04 1.09D-03 4.94D-03 1.61D-02
20 1.10D-09 1.34D-06 6.16D-05 6.88D-04 3.37D-03 1.16D-02
21 3.46D-10 6.21D-07 3.46D-05 4.34D-04 2.29D-03 8.34D-03
22 1.10D-10 2.88D-07 1.95D-05 2.74D-04 1.56D-03 6.00D-03
23 3.46D-11 1.34D-07 1.10D-05 1.73D-04 1.06D-03 4.32D-03
24 1.10D-11 6.21D-08 6.16D-06 1.09D-04 7.25D-04 3.11D-03
25 3.46D-12 2.88D-08 3.46D-06 6.88D-05 4.94D-04 2.24D-03
26 1.10D-12 1.34D-08 1.95D-06 4.34D-05 3.37D-04 1.61D-03
27 3.46D-13 6.21D-09 1.10D-06 2.74D-05 2.29D-04 1.16D-03
28 1.10D-13 2.88D-09 6.16D-07 1.73D-05 1.56D-04 8.34D-04
29 3.46D-14 1.34D-09 3.46D-07 1.09D-05 1.06D-04 6.00D-04
30 1.10D-14 6.21D-10 1.95D-07 6.88D-06 7.25D-05 4.32D-04
31 3.46D-15 2.88D-10 1.10D-07 4.34D-06 4.94D-05 3.11D-04
32 1.10D-15 1.34D-10 6.16D-08 2.74D-06 3.37D-05 2.24D-04
33 3.46D-16 6.21D-11 3.46D-08 1.73D-06 2.29D-05 1.61D-04
34 1.10D-16 2.88D-11 1.95D-08 1.09D-06 1.56D-05 1.16D-04
35 3.46D-17 1.34D-11 1.10D-08 6.88D-07 1.06D-05 8.34D-05
36 1.10D-17 6.21D-12 6.16D-09 4.34D-07 7.25D-06 6.00D-05
37 3.46D-18 2.88D-12 3.46D-09 2.74D-07 4.94D-06 4.32D-05
38 1.10D-18 1.34D-12 1.95D-09 1.73D-07 3.37D-06 3.11D-05
39 3.46D-19 6.21D-13 1.10D-09 1.09D-07 2.29D-06 2.24D-05
40 1.10D-19 2.88D-13 6.16D-10 6.88D-08 1.56D-06 1.61D-05
Table 2: Guide to selecting a value of the discretization paramete
In working with MICROSCOPE, one will rarely be interested in obtaining accurate numerical values of derivatives. Rather, qualitative features will be of primary concern. However, if accuracy is paramount, then the discretization parameter should be picked as small as possible without the round-off term swamping the truncation error. The optimal step size is obtained by equating the round-off term with the truncation error and solving for h. This yields:
(k+2) 1/(k+2)
h = [(1+abs(PHI))*tau*alpha / (gamma *PHI )]
k k
Table 3 lists values of h for the exponential.
Digits (D) Degree of Derivative (k)
1 2 3 4 5 6
1 1.69D+00 2.49D+00 4.50D+00 8.25D+00 1.18D+01 1.88D+01
2 7.83D-01 1.40D+00 2.84D+00 5.62D+00 8.49D+00 1.41D+01
3 3.63D-01 7.87D-01 1.79D+00 3.83D+00 6.11D+00 1.06D+01
4 1.69D-01 4.43D-01 1.13D+00 2.61D+00 4.40D+00 7.93D+00
5 7.83D-02 2.49D-01 7.13D-01 1.78D+00 3.17D+00 5.95D+00
6 3.63D-02 1.40D-01 4.50D-01 1.21D+00 2.28D+00 4.46D+00
7 1.69D-02 7.87D-02 2.84D-01 8.25D-01 1.64D+00 3.34D+00
8 7.83D-03 4.43D-02 1.79D-01 5.62D-01 1.18D+00 2.51D+00
9 3.63D-03 2.49D-02 1.13D-01 3.83D-01 8.49D-01 1.88D+00
10 1.69D-03 1.40D-02 7.13D-02 2.61D-01 6.11D-01 1.41D+00
11 7.83D-04 7.87D-03 4.50D-02 1.78D-01 4.40D-01 1.06D+00
12 3.63D-04 4.43D-03 2.84D-02 1.21D-01 3.17D-01 7.93D-01
13 1.69D-04 2.49D-03 1.79D-02 8.25D-02 2.28D-01 5.95D-01
14 7.83D-05 1.40D-03 1.13D-02 5.62D-02 1.64D-01 4.46D-01
15 3.63D-05 7.87D-04 7.13D-03 3.83D-02 1.18D-01 3.34D-01
16 1.69D-05 4.43D-04 4.50D-03 2.61D-02 8.49D-02 2.51D-01
17 7.83D-06 2.49D-04 2.84D-03 1.78D-02 6.11D-02 1.88D-01
18 3.63D-06 1.40D-04 1.79D-03 1.21D-02 4.40D-02 1.41D-01
19 1.69D-06 7.87D-05 1.13D-03 8.25D-03 3.17D-02 1.06D-01
20 7.83D-07 4.43D-05 7.13D-04 5.62D-03 2.28D-02 7.93D-02
21 3.63D-07 2.49D-05 4.50D-04 3.83D-03 1.64D-02 5.95D-02
22 1.69D-07 1.40D-05 2.84D-04 2.61D-03 1.18D-02 4.46D-02
23 7.83D-08 7.87D-06 1.79D-04 1.78D-03 8.49D-03 3.34D-02
24 3.63D-08 4.43D-06 1.13D-04 1.21D-03 6.11D-03 2.51D-02
25 1.69D-08 2.49D-06 7.13D-05 8.25D-04 4.40D-03 1.88D-02
26 7.83D-09 1.40D-06 4.50D-05 5.62D-04 3.17D-03 1.41D-02
27 3.63D-09 7.87D-07 2.84D-05 3.83D-04 2.28D-03 1.06D-02
28 1.69D-09 4.43D-07 1.79D-05 2.61D-04 1.64D-03 7.93D-03
29 7.83D-10 2.49D-07 1.13D-05 1.78D-04 1.18D-03 5.95D-03
30 3.63D-10 1.40D-07 7.13D-06 1.21D-04 8.49D-04 4.46D-03
31 1.69D-10 7.87D-08 4.50D-06 8.25D-05 6.11D-04 3.34D-03
32 7.83D-11 4.43D-08 2.84D-06 5.62D-05 4.40D-04 2.51D-03
33 3.63D-11 2.49D-08 1.79D-06 3.83D-05 3.17D-04 1.88D-03
34 1.69D-11 1.40D-08 1.13D-06 2.61D-05 2.28D-04 1.41D-03
35 7.83D-12 7.87D-09 7.13D-07 1.78D-05 1.64D-04 1.06D-03
36 3.63D-12 4.43D-09 4.50D-07 1.21D-05 1.18D-04 7.93D-04
37 1.69D-12 2.49D-09 2.84D-07 8.25D-06 8.49D-05 5.95D-04
38 7.83D-13 1.40D-09 1.79D-07 5.62D-06 6.11D-05 4.46D-04
39 3.63D-13 7.87D-10 1.13D-07 3.83D-06 4.40D-05 3.34D-04
40 1.69D-13 4.43D-10 7.13D-08 2.61D-06 3.17D-05 2.51D-04
Table 3: Values of h for maximum accuracy
It turns out that very small jumps (of magnitude eta, say) in some derivative can be detected only by examining a higher order derivative rather than the discontinuous derivative itself. We determine optimal values of h and minimum values of eta by imposing two requirements. First, the delta term should be significantly larger than round-off errors, and, second, the delta term should be significantly larger than the rise of
(k) PHIover half the window. This leads to the equations
m-k -k
eta*delta *h = psi *h *(1+abs(phi))*tau*alpha
m,k 2 k
and
m-k (k+1)
eta*delta *h = psi *h*PHI
m,k 3
where
psi and psi 2 3are safety factors.
We consider two cases. If m = 0 (i.e. there is a step in the trial function) equation (4) is independent of h. Solving for eta yields the minimum value
eta = psi *(1+abs(PHI))*tau*alpha / delta
2 k 0,k
The maximum value of h (above which the delta term will be obscured by the growth in PHI^(k)) is obtained by solving (5) for h:
(k+1) 1/(k+1)
h = [eta*delta / (psi *PHI )]
0,k 3
If m > 0 (i.e. the discontinuity occurs in a derivative) the minimum value of eta can be computed by eliminating h between (4) and (5), and solving for eta. This gives
(k+1) m
eta = [(psi *abs(PHI ))
3
k+1-m 1/(k+1)
*(psi *(1+abs(PHI))*tau*alpha ) )] / delta m,
2 k (
The corresponding optimal value of h (above which the delta term
is too small to be visible in the general growth of PHI^(k) and
below which it is obscured by round-off errors) can be obtained
by eliminating eta between (4) and (5) and solving for h, giving
1/(k+1)
h = (psi *(1+abs(PHI))*alpha *tau / psi )
2 k 3
Notice that this expression is independent of m, the degree of
the derivative in which the jump occurs.
Table 4 lists minimum values of eta (in the last seven columns)
and the corresponding optimal value of h (in the second column)
for the case of the exponential, and for a few values of tau.
Notice how the size of the jump that can be detected is
sometimes decreased by several orders of magnitude if higher
order derivatives are examined.
tau = 1.0D-10
Degree of discontinuous derivative
k optimal h 0 1 2 3 4 5 6
0 5.0D-10 1.5D-09
1 4.5D-05 1.2D-08 1.3D-04
2 1.6D-03 1.2D-08 7.6D-06 4.8D-03
3 1.2D-02 1.8D-08 2.9D-06 4.6D-04 3.7D-02
4 4.8D-02 1.6D-08 1.0D-06 1.1D-04 6.9D-03 1.4D-01
5 1.2D-01 2.0D-08 7.8D-07 3.9D-05 1.6D-03 4.0D-02 3.5D-01
6 2.4D-01 1.9D-08 4.0D-07 1.7D-05 3.7D-04 1.0D-02 8.7D-02 7.2D-01
tau = 1.0D-17
Degree of discontinuous derivative
k optimal h 0 1 2 3 4 5 6
0 5.0D-17 1.5D-16
1 1.4D-08 1.2D-15 4.2D-08
2 7.4D-06 1.2D-15 1.6D-10 2.2D-05
3 2.2D-04 1.8D-15 1.6D-11 1.5D-07 6.6D-04
4 1.9D-03 1.6D-15 2.5D-12 7.0D-09 1.1D-05 5.7D-03
5 7.9D-03 2.0D-15 1.1D-12 8.3D-10 4.9D-07 1.9D-04 2.4D-02
6 2.4D-02 1.9D-15 4.0D-13 1.7D-10 3.7D-08 1.0D-05 8.7D-04 7.2D-02
tau = 1.0D-24
Degree of discontinuous derivative
k optimal h 0 1 2 3 4 5 6
0 5.0D-24 1.5D-23
1 4.5D-12 1.2D-22 1.3D-11
2 3.4D-08 1.2D-22 3.5D-15 1.0D-07
3 3.9D-06 1.8D-22 9.1D-17 4.6D-11 1.2D-05
4 7.6D-05 1.6D-22 6.3D-18 4.4D-13 1.7D-08 2.3D-04
5 5.4D-04 2.0D-22 1.7D-18 1.8D-14 1.6D-10 8.7D-07 1.6D-03
6 2.4D-03 1.9D-22 4.0D-19 1.7D-15 3.7D-12 1.0D-08 8.7D-06 7.2D-03
tau = 1.0D-31
Degree of discontinuous derivative
k optimal h 0 1 2 3 4 5 6
0 5.0D-31 1.5D-30
1 1.4D-15 1.2D-29 4.2D-15
2 1.6D-10 1.2D-29 7.6D-20 4.8D-10
3 7.0D-08 1.8D-29 5.1D-22 1.5D-14 2.1D-07
4 3.0D-06 1.6D-29 1.6D-23 2.8D-17 2.8D-11 9.1D-06
5 3.7D-05 2.0D-29 2.5D-24 3.9D-19 4.9D-14 4.0D-09 1.1D-04
6 2.4D-04 1.9D-29 4.0D-25 1.7D-20 3.7D-16 1.0D-11 8.7D-08 7.2D-04
Table 4: The superiority of testing Delta Functions and their derivatives
This appendix contains (1) a brief discussion of the organizational framework upon which MICROSCOPE is built, (2) a detailed explanation of the sampling scheme and logic which MICROSCOPE uses to calculate the tangential derivatives, and (3) a brief description of the common blocks used.
MICROSCOPE has evolved over a period of several years, progressing through several preliminary versions. The core of the present version was designed and written by Bill Harris during the summer of 1983. MICROSCOPE then grew, under the hands of Peter Alfeld, to a substantially larger package which is the current version. This explains lack of uniformity in the commenting of the source code. It also accounts for the inconsistency in the notation used within the program versus that employed in this manual. However, only very serious users interested in understanding and possibly manipulating the code (for whome this appendix is written) should suffer from such idiosyncrasies of the package.
MICROSCOPE follows a common and well established format for interactive programs in FORTRAN. The user must provide a short main program, which assigns FORTRAN Input/Output device numbers and screen parameters, and then calls a master routine (MCRSCP). MCRSCP is a driver program from which all other routines are called and to which execution returns after each command has been completed. This is accomplished by (i) prompting for a new command, (ii) reading and recognizing the new command, (iii) using a computed GO TO to branch to the appropriate routine, (iv) a bit of logic and computation (i.e. sampling the function and calculating the derivatives), (v) displaying the results in some graphic form, and (vi) updating an options table and returning to step (i).
To help understanding the following flowchart several ideas and their representative variables need to be introduced. The logical variable LGO indicates whether the screen is active (LGO = .TRUE) or inactive (LGO = .FALSE.). It is set and reset at a number of places in the package. The only place where the value of LGO has an effect on the program flow is in the routine MCRSCP, between statements 7900 and 8300.
There are 4 basic types of commands. These correspond to the integer variable ITYPE which is set appropriately by each command. Commands of ITYPE = 1 imply that the function must be sampled from scratch. Since other changes may be desired before computation begins, these commands deactivate the screen and this is accomplished by setting LGO = .FALSE.. Commands which are of ITYPE = 2 require some computation. Two special ITYPE = 2 commands are FORCE and GO. These cause LGO to be set to .TRUE. and are the only commands that activate an inactive screen. The remaining commands of ITYPE = 2 only make modifications of the current setup and usually only require small amounts of additional computation. When LGO = .TRUE. any ITYPE = 2 command is executed before the next command is prompted for. ITYPE = 3 commands do not require computation but cause the alphanumerical display to be altered. Finally, commands of ITYPE = 4 do not alter the display and do not require computation.
The other major logical control variable is LSCRN. Its main use is in screen-editing versions of MICROSCOPE. LSCRN = .FALSE. means that the alphanumerical display has been overwritten by text and must be regenerated. LSCRN = .TRUE. means that the display is intact from the previous command and only selected pixels need to be updated. LSCRN is used by the display routine SGRAPH in the file SCREEN .
Computation time for evaluating display functions can be sizable. The program logic and sampling scheme were designed to utilize previous computations and eliminate unnecessary recomputations to the extent possible. For example, shifts and magnification changes often require only a few additional evaluations.
The flowchart that follows gives an overview of the logical flow in the MCRSCP driver routine. It concentrates on the portion between statement 1500 and 8300 where the logic is somewhat involved.
Epitome Flowchart
!
! from Main Program
!
V
* output message identifying MICROSCOPE
* check screen parameters
* copy channel numbers, screen parameters, and set flags
* set LGO = .FALSE.
* compute round-off characteristics (NUMDIG)
* set default values for some of the variables contained in
the various common blocks
* input some of the help documentation as well as news (INHELP)
1500 * prompt/input/parse command (DIALOG)
* reset error flags
* branch to the appropriate subroutine (using a computed GO TO)
This corresponds to statements 1600 through 7800. Many
of these commands alter variables, with subscript 2, in
the common block OPTION. These will then be used by
later routines.
7900 * if command is of ITYPE = 1 (i.e. requires complete recomputation)
then set: LGO = .FALSE.
* if ITYPE = 2 and LGO = .TRUE.
then calculate values for NUM and DEN corresponding to
sampling changes. (See documentation in part 2
of this appendix) (SGAMMA)
* check for errors and output appropriate messages if any occurred
* if error = .TRUE.
<------------then go to 1500
* if LGO = .FALSE. and the Refresh Screen (RSCREEN) command is given
then display the old graphical information (SGRAPH)
* if LGO = .FALSE.
<------------then go to 1500
* branch to the appropriate section of the program depending on
the value of ITYPE (using a computed GO TO)
8000 * (if ITYPE = 2)
* check to see which aspects of the sampling have been altered
by the command and set the appropriate logical flags (CHKCMP)
* call the routines flagged by the CHKCMP routine.
(the actual sampling and computation of derivatives occurs
only in the routine SCMUPD)
8100 * (if ITYPE = 3)
* set up the graphics data (i.e. update plot array) (SGRUPD)
* graph the data (SGRAPH)
* if logging is required
then output a copy of the graphical display in scroll mode
* copy current set of parameter values from new to old. (OPCOPY)
<-------* go to 1500
In the remainder of this appendix an explanation of the sampling scheme and computational logic is given. The variable h, as used in the source code, is not the same as the variable labeled as h on output. To reduce confusion we will use the following notation in this appendix (but nowhere else in this manual):
Variable in Variable in Variable in Output
Source Code Appendix and elsewhere in
the manual
h ss (step size) s
sw sw (stencil width) h
The output names correspond to the common usage while the names used in the program and appendix have the meanings defined below.
The basic idea of the sampling scheme consists of representing the relevant part of the line of investigation in two linear arrays, one corresponding to points in the domain, and one containing corresponding values of the display function (if it has been evaluated). Any changes in the stencil width and step size, as well as small shifts of the point of examination are then implemented by shuffling the information in these two arrays. The size of the arrays has been chosen to accomodate fairly large changes in parameter values with little computation. The key variables and arrays used in the program are as follows:
The sampling points are spaced unevenly along the line of investigation, but evenly in index space. In the following figures, asterisks denote the points at which the display function needs to be evaluated. A hyphen marks a point for which space for a function value has been provided in the arrays XS and FS, but where no evaluation is necessary. Consider a point p at which we wish to evaluate the display function and its derivatives up to order 6. The stencil is centered about p.
p
-.......-...-...-.......*.......-...-...-.......- f(p)
.
. (1)
*.......-...-...-.......-.......-...-...-.......* f (p)
.
. (2)
*.......-...-...-.......*.......-...-...-.......* f (p)
.
. (3)
*.......-...*...-.......-.......-...*...-.......* f (p)
.
. (4)
*.......-...*...-.......*.......-...*...-.......* f (p)
.
. (5)
*.......*...-...*.......-.......*...-...*.......* f (p)
.
. (6)
*.......*...-...*.......*.......*...-...*.......* f (p)
sw sw sw sw sw sw sw sw
- -- - -- - -- - -- 0 -- -- -- --
2 3 4 6 6 4 3 2
<---------------------- sw --------------------->
(h on output)
The figure describes the actual sampling structure along the line of investigation. Space has been provided only for points at which the display function may have to be evaluated. Providing space for sampling at increments of sw/12 would be simpler but wasteful since at most only 9 points (not 13) are needed in each stencil.
Next, consider the problem of extending the sampling over an entire line segment. The particular extension used in MICROSCOPE was chosen so that changing the step size and stencil width by powers of two would require minimal computation. Other factors are allowed but require greater computation.
Considering a portion of the line segment containing a sequence of three plotted points: p1,p2, and p3, the spacing (for the ratio sw/ss = 2/1, nw = 4, and assuming all tangential derivatives are to be plotted) looks like:
..*.....*..*..*.....*.....*..*..*.....*.....*..*..*.....*.....*..*..*.....*..
. p1 p2 p3 .
. . . . .
*.....*..*..*.....*.....*..*..*.....* (stencil for p1). .
-4 -3 -2 -1 0 1 2 3 4 . .
. . . .
*.....*..*..*.....*.....*..*..*.....* (stencil for p2).
-4 -3 -2 -1 0 1 2 3 4 .
. . . .
.(stencil for p3) *.....*..*..*.....*.....*..*..*.....*
-4 -3 -2 -1 0 1 2 3 4
<---------------- sw --------------->
<------ ss ------> <------ ss ------>
Note how the stencils overlap, allowing the use of one value of the display function for more than one (here two or three) values of a derivative. If the stencil width is doubled while the step size is kept constant (corresponding to an increase in the window width by a factor 2) we obtain:
..*.....*..*..*.....*.....*..*..*.....*.....*..*..*.....*.....*..*..*.....*..
. p1 p2 p3 .
..*.....*..-..-.....*.....-..-..*.....*.....*.....-.....* (stencil for p1).
-4 -3 0 3 4 5 8 .
. . . . .
*........-..*.....*.....*..-..-.....*.....-.....*.....*.....*...........*
-8 -5 -4 -3 0 3 4 5 8
. . . .
.(stencil for p3) *.....-..-..*.....*.....*...........*.....-.....*.....*..
-8 -5 -4 -3 0 3 4
. . . . .
<---------------------------------- sw --------------------------------->
<------ ss ------> <------ ss ------>
The above operation requires additional evaluation of the display function only at the boundary of the interval of investigation. The reverse process of halving the window width would require no additional evaluations at all.
So that all derivatives, up to and including the sixth, can be evaluated, the minimal number of points between plotted points (such as p1,p2,p3) is 4. Furthermore, the distance is always a multiple of 4. Similarly, the minimal number of points needed in a stencil is 9 (i.e. -4,...,0,...,4), and it is always of the form 1+(an integer multiple of 8), although evaluation takes place at at most 9 points. These features are illustrated in the following progression of increasing stencils:
..*.....*..*..*.....*.....*..*..*.....*.....*..*..*.....*.....*..*..*.....*.. . . p2 . . . *.....*..*..*.....*.....*..*..*.....* . . -4 -3 -2 -1 0 1 2 3 4 . . . . . . *.....-..-..*.....*.....*..-..-.....*.....-..-..*.....*.....*.....-.....* -8 -5 -4 -3 0 3 4 5 8 . . . . . ..*.....-..*..-.....*.....-..-..-.....*.....-..-..-.....*.....-..*..-.....*.. -8 -6 -4 0 4 6 8
These 3 stencils span index sets of 9,17, and 25 points, respectively (although the third set would not fit on this page). The variable nw for these is 4,8, and 12, respectively.
The relation between the quantities ss and ROPSTS(2) is given by
nh
ss = -- * ROPSTS(2)
8
The program strives to maintain a value of nh = 8. For a more detailed discussion see the description of the routine SGAMMA below.
The relation between the point of examination, P, a point of evaluation, PS, the distance t of ps from P, the direction of examination, d, and the quantities in MICROSCOPE, is given by:
.
.
.
XS(ICENTR-2) = -ROPSTS(2)/4
XS(ICENTR-1) = -ROPSTS(2)/6
XS(ICENTR) = 0
XS(ICENTR+1) = ROPSTS(2)/6
XS(ICENTR+2) = ROPSTS(2)/4
XS(ICENTR+3) = ROPSTS(2)/3
XS(ICENTR+4) = ROPSTS(2)/2
XS(ICENTR+5) = 2*ROPSTS(2)/3
XS(ICENTR+6) = 3*ROPSTS(2)/4
XS(ICENTR+7) = 5*ROPSTS(2)/6
XS(ICENTR+8) = ROPSTS(2)
.
.
.
and
P = ROPPNT(*,2)
d = ROPUDI(*,2)
t = XS(i)
PS = P + t * d
FS(i) = f(PS)
LDEF(I) = .TRUE. <===> FS(I) has been evaluated
To explain the calculation of the derivatives, consider a point p at which we wish to evaluate the function and its derivatives. To evaluate all derivatives up to the sixth we need 9 points. We number them as follows:
1 2 3 4 5 6 7 8 9
*.......*...*...*.......*.......*...*...*.......*
. . . . p . . . .
. . . . . . . . .
sw sw sw sw 0 sw sw sw sw
- -- - -- - -- - -- -- -- -- --
2 3 4 6 6 4 3 2
Instead of using the formulas given in section 1.2 directly, the following equivalent but more efficient formulas are used:
g0 = f5
g1 = f1 + f9
g2 = - ( f1 - f9 )
g3 = 4 * ( f3 + f7 )
g4 = 4 * ( f3 - f7 )
g5 = 4 * ( f2 - f8 ) + 5 * ( f6 - f4 )
g6 = - 6 * ( f2 + f8 ) + 15 * ( f4 + f6 )
(0)
Q = g0
i
(1)
Q = (1/sw) * g2
i
(2) 2
Q = (2/sw) * ( g1 - 2*g0 )
i
(3) 3
Q = (4/sw) /4 * ( 2*g2 + g4 )
i
(4) 4
Q = (4/sw) * ( g1 - g3 + 6*g0 )
i
(5) 5
Q = (6/sw) /2 * ( g2 + g5 )
i
(6) 6
Q = (6/sw) * ( g1 + g6 - 20*g0 )
i
The sequence in which MCRSCP calls the main routines is:
SGAMMA - calculates gamma (see following discussion)
CHKCMP - checks to see what changes in computations are needed
SSMPUD - shift and copy
SCMUPD - compute new values
When a step size/stencil width changing command is given, the program will attempt to change nh and nw without having to change ROPSTS(2) (and hence copy or throw out data). Otherwise, MCRSCP attempts to choose a change in scaling (ROPSTS(2)) which saves as much data as possible. This is the purpose of the routine SGAMMA which we now discuss in some detail. First we introduce additional notation.
Primes will denote the new values for nh,nw,ss,and s. For example: h' denotes the step size after changes have been made.
Define: gamma = s/s' = num/den > 0
Constraints: The commands set the values of m and v. There are 5 constraints which are then used to choose the value of gamma. They are not all independent of each other.
ss = (nh/8) * s and ss' = (nh'/8) * s'
Also:
ss' = m * h
Hence:
ss' = (nh'/8)*s' = m * h = m * (nh/8) * s
===> (nh'/8) * s' = m * (nh/8) * s
===> nh' = m * nh * s/s' = m * nh * gamma
Similarly:
sw = (nw/4) * s, sw' = (nw'/4) * s', and sw' = v * sw
===> sw' = (nw'/4) * s' = v * sw = v * (nw/4) * s
===> nw' = v * nw * s/s' = v * nw * gamma
In summary:
nh' = m * nh * gamma
nw' = v * nw * gamma
Referring to the conditions numbered above, divide (1) by m * nh
and (2) by v * nw to get:
4 IXSIZE 1
-------- <= gamma <= ------- * --------
m * nh width m * nh
and
4 IXSIZE 1
-------- <= gamma <= ------- * --------
v * nw 2 v * nw
Let
4 4
a = max( ------- , ------- )
m * nh v * nw
and
IXSIZE 1 IXSIZE 1
b = min( -------- * ------- , -------- * -------- )
width m * nh 2 v * nw
From above, gamma must satisfy:
a <= gamma <= b
Furthermore:
nh' = m * nh * gamma and mod(nh',4) = 0
===> (*) mod( m * nh * gamma , 4 ) = 0
nw' = v * nw * gamma and mod( v * nw * gamma , 4 ) = 0
===> (**) mod( v * nw * gamma , 4 ) = 0
The routine SGAMMA uses the above two conditions (*) and (**) to choose acceptable values for gamma (and hence num and den since gamma = num/den).
* check to see if the changes m and v can be accomodated by changing
nh and nw, but not ss (h in code).
* if no changes in ss are needed
<---------- then num = den = 1 ==> gamma = 1 and return
* calculate a and b
* if m = v
<---------- then set gamma = m, ROPSTS(2) = m * ROPSTS(2), and return
* calculate rho
m' = rho * m (m' and v' are positive integers)
v' = rho * v
* determine "all" possible numerators and denominators for gamma
(up tp 32) which satisfy:
a <= num/den <= b
mod( m * nh * gamma , 4 ) = 0
mod( v * nw * gamma , 4 ) = 0
* rank the candidate values of gamma according to 4 criteria
(1) nh' * width + 2 * nw' <= IXSIZE
(2) 8 * width <= nh' * width + 2 * nw' <= ICENTR
(this allows for halving and doubling without rescaling)
(3) gamma a power of 2 (positive or negative)
(4) nh' > 4 and nw' > 4
* choose the combination of num and den which has the highest cor-
responding rank and set the variables accordingly.
<--* return
Here we list the COMMON blocks used in MICROSCOPE. This information may be useful in modifying the code and in accessing the data structures of the package via the command USER (see section 4.15). For each common block we list in parentheses the number of words (i.e. 2 words for each DOUBLE precision variable and 1 word for each REAL, INTEGER, or LOGICAL variable) occupied by the block. The blocks are listed in the same sequence as they appear in the routine MCRSCP.
CB (16) controls cross differentiation.
CBDSW : width of differentiation stencil
CBD(3) : normalized cross direction
CBDU(3) : unnormalized cross direction
ICBD : order of cross differentiation
LCBD : .TRUE. if cross differentiation is on, .FALSE.
otherwise
IO (5) contains the input/output device numbers.
RECORD : Record device number. (Output data file for recording
screen images and comments)
GRAPHD : Graphics device number. (terminal or graphics display)
HELPD : Help file device number.
RSTRTD : Restart file device number. (unformatted output data
file for use in restarting at a later date)
INPUTD : Input device number. (usually the terminal)
LOG (1) controls logging
LCHN : Logging device number (= 0 <===> logiing is off
LOGCOM (5391) contains logical arrays used primarily by SCMUPD
to keep track of what has already been computed.
LDF(7) : Denotes whether the (i-1)'th derivative has been calcu-
lated. LDF(i) = .FALSE. means the (i-1)'th derivative
has not been calculated. LDF(i) = .TRUE. means it has.
LPLT(7) : Denotes whether the (i-1)'th derivative is currently
being displayed. LPLT(i) = .FALSE. means the (i-1)'th
derivative is not being displayed and LPLT(i) = .TRUE.
means it is.
LDEF(5377) : Denotes whether or not the interpolation function has
been calculated at the coordinates corresponding to
XS(i). (i.e. whether FS(i) has been assigned)
LDEF(i) = .FALSE. means the i'th value has not been
sampled and LDEF(i) = .TRUE. means it has.
FUNCOM (23,426) contains arrays pertaining to the values of the
interpolation function at its sampled points.
XS(5377) : The distance array, always centered about 0. Its defini-
tion and usage is given in the second part of this
appendix.
FS(5377) : The array containing the values of the interpolation
function corresponding to the points specified by the
XS() array. See the second part of this appendix for
discussion.
DF(135,7) : The finite differences. DF(i,j) corresponds to the i'th
plotted point and (j-1)'th derivative. These are the
p1,p2,p3's of the discussion in the second part of
this appendix.
DFMNMX(2,7): These are the minimums/maximums of the finite differences
in the DF() array.
ERRCOM (2) used to handle error checking and reporting.
When errors are detected in any of the sub routines
or computation routines, the errors are flagged and
then reported by the CHKERR routine.
ERRCOD : This is an integer code number for the type of error
detected.
ERR : This is the logical flag which denotes whether or not
an error has occurred. ERR = .FALSE. means no error
has occurred and ERR = .TRUE. means an error has
occurred.
PLTCOM (16342) contains display versions of the derivatives.
IPLOT(135,7): IPLOT(i,j) is the "screen coordinate" (integer) of the
(j-1)'th derivative in the i'th column.
ISCRN1(135,57): The previous screen image. ISCRN1(i,j) contains the
Hollerith representation of the character to be
displayed at the (i,j)'th pixel on the screen.
(where indexing starts at the upper left corner
of the screen). i is the horizontal component and
j is the vertical component.
ISCRN2(135,57): The current screen image. This has the same description
as ISCRN1(,) above.
SCALE(7) : The scale factors used to individually scale, up or down,
the finite difference "curves" (in the DF(,) array)
so that the curve fills up the entire screen in the
sense that the maximum point lies on the top row and
minimum point lies on the bottom row of the screen.
LOADIO (3) contains channel numbers of the loading devices.
PLOAD : point loading device
DLOAD : direction loading device
CLOAD : cross direction loading device
SCREEN (7) contains the values of parameters describing the
device on which most interaction occurs (i.e. the
screen of the terminal in most cases)
OUTPUT : The device number of the graphics device (usually the
terminal)
LINES : The number of usable lines on the screen. Note that
this may have to be one less than the number of
displayable lines.
WIDTH : The number of columns to be displayed. It is best to
choose an odd number for this so that the center
can be plotted. It must be less than or equal to
135. 75 is a common value.
ILP : The number of lines in the curve plotting portion of the
display. It must be less than or equal to 57.
IPRMPT : The number of lines above the screen bottom at which
the prompting commands are given. It must be at
least 2.
IDSPLA : The number of lines above the screen bottom at which the
numerical data region begins. It must be greater than
IPRMPT by 6 or more.
LSCRN : The logical variable denoting whether the curves are
currently being displayed or they have been over-
written by text. This is to allow the program to
take advantage of terminals which allow screen
editing. LSCRN = .FALSE. means that the screen has
been overwritten since the previous command and
LSCRN = .FALSE. means that the screen is intact from
the previous command.
OPTION (122) contains most of the variables, other than the
sampling arrays in FUNCOM, which may change due to
commands. There are 2 copies of each variable.
(1) is the previous value and (2) is the current
value. If mistakes are made in the inactive (deactiv)
mode (i.e. LGO = .FALSE.) then the UNDO command will
erase the changes and copy the old set (1) to (2).
ROPDI(3,2) : The vector specifying the direction of investigation.
Note that for all vector valued quantities, the
dimension is specified by IOPDM() below.
ROPPNT(3,2): The vector specifying the center of the line of
investigation. If the shift is 0 (i.e. IOPSH(2) = 0)
then this is the point which is at the center of the
screen.
ROPDR1(3,2): Currently unused vector.
ROPDR2(3,2): Currently unused vector.
ROPSTS(2) : The nominal spacing between the points at which the
finite differences are calculated. These will be
the plotted points. "Nominal" refers to the base
step size when nh(2) = 8, as discussed in the second
part of this appendix. As noted earlier, this is
not the h which appears on output.
ROPSTW(2) : The stencil width (sw). This corresponds to the h which
appears on output. It is discussed in the second
part of this appendix.
ROPUDI(3,2): The normalized (unit) vector specifying the direction
of investigation. This is the normalized version
of the ROPDI() vector.
NH(2) : The step size in index space. See the discussion in the
second part of this appendix.
NW(2) : Half of the stencil width in index space. See the
discussion in the second part of this appendix.
ILEFT(2) : The index of the left-most element of the FS() array
which has been assigned a value.
IRIGHT(2) : The index of the right-most element of the FS() array
which has been assigned a value.
IOPDM(2) : The dimension of the domain space. It can be 1,2 or 3.
IOPSH(2) : The shift, in integer increments, of the plotted points
from the center defined by ROPSTS(*,2).
IOPSS(2) : The scaling factor for the step size. It is set by
the commands.
IOPWN(2) : The scaling factor for the stencil width. This
multiplies/divides the current window, as it appears
on the screen. It is set by the commands.
OPDC(2) : The logical flag determining whether or not the center
of the screen is to be marked.
OPDS(2) : The logical flag determining whether or not a scale is
to be displayed along the horizontal axis.
OPDX(2) : The logical flag determining whether or not a line of
dashes will be displayed to represent the x-axis.
OPCMP(2) : The logical flag denoting whether or not any computations
have yet been done for the current set of ROPPNT and
ROPDI.
OPSPL(2) : Currently not in use.
OPDF1(7,2) : The logical array specifying which derivatives are to
be displayed. It is set by the DGRAPH and EGRAPH
commands.
OPDF2(7,2) : The logical array specifying which derivatives are to
to be accentuated as they are displayed. It is set
by the ACCENT command.
NRMLZE (1) controls normalization of derivatives
NORMAL : = .TRUE. <===> normalization is on
ROOM (2) controls whether one line of the graphical
display is transferred to the numerical display.
This is only necessary if all six tangential
derivatives of a cross derivative are plotted.
ILPUSR : User specified space for the graphical display
IDSUSR : User specified space for the numerical display
The corresponding variables ILP and IDSPLA in
the common block SCREEN are modified as
necessary.
HELPER (7330) contains help information
HELP : Help device channel number
JHELP1 : Length of the HSUMRY information (in lines)
JHELP2 : Length of the HELP information (in lines)
JHELP3 : Number of available commands
JHELP(99,2): (I,1) line in which detailed help for I-th
command begins
(I,2) length of detailed help
IHELP(72,99): (*,I) prompt text for I-th command
USER (5) see appendix III
PLOWN (296) contains current PLOT options
ITITLE(72) : Title
IBOTTM(72) : Legend
NUMBR : page number
BL(20) : left bottom label
BR(20) : right bottom label
FRAME : = .TRUE. <===> frame is drawn
COLOR : = .TRUE. <===> plot is in color
NUMRCL : = .TRUE. <===> numerical display is drawn
LMARK : = .TRUE. <===> graphs are marked
LBLS : = .TRUE. <===> graphs are labeled
DATE : = .TRUE. <===> date is drawn
TIME : = .TRUE. <===> time is drawn
FOOWN (1) controls revaluation after the command FORCE has
been issued.
LFO : = .TRUE. <===> reevaluate
The examples in this manual can be reproduced, and experience with MICROSCOPE can be gained, by using the test routines (a special choice of the trial function F and the user intervention routine SUBUSR, plus a few support routines) supplied with the package. The test code is contained in the files USRLC and USRUC (see chapter 5).
Several different trial functions are available for examination, and rounding to any number of digits up to the number supported on your local installation can be simulated. All options are selected using the command USER, which puts you
into user mode. You can change options similarly as in the PLOT mode until you are satisfied, and then return to the command mode. Options are chosen by passing numbers to the program according to the prompts given. A current list of options is maintained on the screen.
The trial function is selected by giving its reference number. The parameter eta is 1 by default and can be changed by typing -3. By default rounding to 10 digits is simulated. Rounding can be changed by typing -1 and giving the desired number of significant digits. It can be turned off by giving 0 as the number of digits. The following trial functions are available:
x
1: F(x) = s(0,eta,x) or F(x) = s(0,eta,x) + e
x
2: F(x) = s(1,eta,x) or F(x) = s(1,eta,x) + e
x
3: F(x) = s(2,eta,x) or F(x) = s(2,eta,x) + e
x
4: F(x) = s(3,eta,x) or F(x) = s(3,eta,x) + e
x
5: F(x) = s(4,eta,x) or F(x) = s(4,eta,x) + e
x
6: F(x) = s(5,eta,x) or F(x) = s(5,eta,x) + e
The addition of the exponential in the above functions is off by
default, and is toggled on or off by typing -2.
eta*x
7: F(x) = e
2
8: F(x,y) = eta*abs(x)*x + (1-eta)*abs(y)*y
2
9: F((x,y) = x *abs(x)*y*abs(y)
x
10: F(x) = the spline interpolant to e used in section 2
x
11: F(x) = e - function 10
(i.e. F(x) is the error in the spline 10)
12: F(x) = 0
13: F(x) = eta*x
e
14: F(x) = x ta
e
15: F(x) = x ta
The user mode is left by typing zero (or just RETURN).
The commands STORE and RESTART are effective on the test parameters which are stored in the COMMON block USER. The parameters have the following meanings:
ETA: DOUBLE PRECISION variable containing the value of eta.
IROUND: Number of digits to which function values are rounded (IROUND = 0 <===> no rounding is taking place).
N: reference number of the trial function.
ADD: LOGICAL variable that indicates whether an exponential term is to be added to the trial functions defined by N = 1,...,6.
The following table gives the USER defaults and the parameter settings necessary to generate the figures in this manual. N is the reference number of the trial function, D is the number of digits rounded to. The entry "--- " means "does not apply". The default settings are effective on the first USER and after giving the command SETDF. Other parameters necessary for reproducing the examples in this manual (point, direction, dimension, interval etc.) can be read off the alphanumerical display in each figure.
Figure N D ETA exponential added
Default 1 10 1 NO
1 9 10 --- ---
2-6 11 10 --- ---
7-8 10 10 --- ---
9 12 10 --- ---
10 13 10 2D-9 ---
11-12 10 10 --- ---
13-18 3 10 -5D-3 YES
19 1 10 -1 NO
20 2 10 -1 NO
21 3 10 -1 NO
22 4 10 -1 NO
23 5 10 -1 NO
24 6 10 -1 NO
25-26 7 10 1 ---
27 9 10 --- ---
28 14 10 1.5 ---
29 15 10 1.5 ---
Table of test parameters
Because of the random nature of round-off errors, somewhat different displays may be obtained even with the above parameter settings.
Some of the above test functions so far have not been used in this manual. They illustrate more subtle phenomena that will not normally occur in multivariate interpolation and approximation. In a typical application, the trial function will have a finite degree of differentiability across the front, and be infinitely often differentiable in any direction contained in the front. Function 8 illustrates a situation where the degree of differentiability is finite in any direction, but is dependent upon the direction of investigation. The function is once differentiable in the y direction and twice differentiable in the x direction. The relative sizes of the jumps in the appropriate derivatives can be controlled by the parameter eta.
The functions 14 and 15 can be used to illustrate discontinuities other than step discontinuities. The functions differ in the way they are defined for negative values of x, for which the value of function 14 is negative, and that of function 15 positive, if the integer part of abs(eta) is even, and of the opposite sign if it is odd.
For both functions 14 and 15, the value at x = 0 is 0 if eta is negative. This prevents a floating point overflow and a meaningless display.
The following two figures illustrate functions 14 and 15 for eta = 1.5 (i.e. a function that is 1.5 times differentiable). The first derivative is continuous, but the second derivative exhibits a pole. The illustration of these features is enhanced by choosing the window width to be less than 1.
.. 2 1111111111
.. I 11111111 ..
.. I 111111 ..
. I 111111 .
.. I 111 ..
.. I 111 ..
.. I1 ..
.. 1 ..
... 1I ...
.. 111 I ..
.. 1111 I ..
..11111 2I2 ...
1111111 ... 2 I 2 ...
11111111 ...2222 I 2222...
2222222222222222222222222222222 ........... 2222222222222222222222222222222
===========================================================================
Point = 0.000000D+00 s = 4.0000D-04
Direction = 1.000000D+00 h = 2.0000D-04
F0 ( 2.14D-11, 1.80D-03) F1 (-1.82D-01, 1.82D-01) F2 ( 6.17D+00, 1.41D+02)
I/O: 25 6 6 1 2 3 27 28 29 30 NRML on current CALLS = 2901
Figure 28 A 3/2 times differentiable function
11111 I2 11111
1111 I 1111.
1111 I 2 1111..
1111 I 22 1111..
111 I 2222 111..
111 I 22222222222222222
111 I ......111 222222222222
11 ......2..... 11
222222222222222 ...11. I 11
22222222222222 11 I 11
.... 22221 I 11
.... 22 I 1
.... 2 I 1
... 1I1
.... 21
===========================================================================
Point = 0.000000D+00 s = 4.0000D-04
Direction = 1.000000D+00 h = 2.0000D-04
F0 (-1.80D-03, 1.80D-03) F1 ( 1.41D-02, 1.82D-01) F2 (-3.81D+01, 3.81D+01)
I/O: 25 6 6 1 2 3 27 28 29 30 NRML on current CALLS = 3052
Figure 29 Another 3/2 times differentiable function
Figure 29 shows a cusp in the first derivative whereas in figure 28 the first derivative passes vertically through the origin. Both derivatives are continuous at the point of examination and have a vertical tangent there. The fact that the second derivative has a pole can be verified by observing that the extreme function value increases as the discretization parameter is decreased. For a cusp, the extreme function value approaches a limit as h decreases. This phenomenon can be utilized to distinguish between a cusp and a pole.
It is of course straightforward to program other or additional examples for a trial function by modifying the routines in the test package appropriately.
MICROSCOPE was developed over several years and has evolved through several preliminary versions. This process was shaped by the work of the Computer Aided Geometric Design Group at the Mathematics Department of the University of Utah and has benefited from the interactions with Paul Arner, R.E. Barnhill, Nelson Beebe, Gerald Farin, Tom Jensen, Chip Petersen, Bruce Piper, Tracy Whelan, Andrew Worsey, and others.
This work was supported partly by the Department of Energy under Contract No. DE-AC02-82ER12046 A000 to the University of Utah, by the United States Army under Contract No. DAAG29-80C-0041, by two grants from the University of Utah Research Committee, and by a sabbatical leave of the first author from the University of Utah. A large part of the program development, and this documentation, were carried out at the Mathematics Research Center at the University of Wisconsin. Fred Sauer of MRC has been most helpful with all computer related problems.
Fine print, your comments, more links, Peter Alfeld, PA1MI
[16-Aug-1996]