Interpolation


Assignment :

Start Homework 11 and Study for the Exam (Wednesday, April 2, 8:15, 102 Forum)

New Fortran:

EXTERNAL and INTRINSIC Statements, SAVE

Digression to EXTERNAL and INTRINSIC type statements

There are a few statements worth knowing that I skipped earlier to focus on more important features. EXTERNAL and INTRINSIC are not frequently used, but every now and again you are glad that you have them available. The example plot2.f uses both EXTERNAL and INTRINSIC type statements. If you look at the program you will see that I actually pass the name of a function through the argument list of the subroutine "gnuplot". This is consistent with the way argument lists work. You should recall that when a variable name is placed in the calling argument list, what happens is that the memory location of the variable is passed to the subprogram. When a function name appears in the argument list the starting location of the that function is passed. As long as the subprogram knows that particular argument is a function, it knows where to branch in memory to start executing the function's instructions. One way that you can clearly identify a given dummy argument as a function is to include the dummy argument's name in an EXTERNAL type statement within the subprogram. As a simple example:

      subroutine dotrig(value, angle, trigfunc)
      implicit none
      real value, angle, trigfunc
      external trigfunc
      value = trigfunc(angle)
      return
      end
However, the "external trigfunc" is not mandatory (except on a few impoverished compilers that I have seen). The fact that "trigfunc" is not declared and array means that "trigfunc(angle)" is interpreted as a the use of a function.

Another rule is that the program unit doing the call with a function in the argument list must clearly declare that to be a function (parentheses are not next to the name when it is used in the argument list). If the function to be passed is a Fortran Intrinsic Function, it should be declared with an INTRINSIC type statement. If the function to passed is another portion of the program itself it should be declared with an EXTERNAL type statement.

      program test
      external sumtrig
      intrinsic sin, cos
      real angle, value1, value2, value3
      data angle/ 1.0/
      call dotrig (value1, angle, sin)
      call dotrig (value2, angle, cos)
      call dotrig( value3, angle, sumtrig)
      stop
      end
      function sumtrig(ang)
      real ang, sumtrig
      sumtrig = cos(ang) + sin(ang)
      return
      end
The need for both an EXTERNAL and INTRINSIC statement is a result of a second use of the EXTERNAL statement. I can bypass any Fortran intrinsic function, by writing my own function with that name, and in every program unit using the function with that name, including an EXTERNAL statement with the name of the replaced function. For example, I could write my own version of the SIN function, and then include a line "external sin" in any program unit (main program, subroutine, or function) that uses "sin". Always remember that EXTERNAL and INTRINSIC are non-executable statements and belong in the top section of the routine with other type statements, before any executable statement.


Interpolation

Remember the interpolation shown in the IF samples inter1.f and inter2.f. This is the same sort of thing but now we are using interpolation to fill in connecting points of data, not simply to smooth the transition between two different regions.

Why bother? One good reason is the use of an equation of state. For an ideal gas you can write a very simple equation of state, but for other gases or fluids it is not so simple. Take steam as an example. Behavior, particularly near the saturation line, is fairly complex. You can use data as a starting point for some very complex curve fits, or use some very complicated calculations. Either approach may require too much computer time. One alternative is to store values of pressure at carefully selected temperatures and densities, then interpolate to get pressures at other values of temperature and density. I'm getting ahead of myself here since we will need a doubly dimensioned array to pull off this particular example. However, one purpose of demonstrating interpolation here is as a lead-in to multiply-dimensioned arrays.

Interpolation is an exercise in curve (or surface) fitting. For linear interpolation, you are fitting a straight line to two points. Look at two points with x, y value pairs (x1, y1 ) and (x2, y2). What is the equation giving y as a function of x for the line passing through these points? At times this simple approach will result in unacceptably large errors (fit versus "truth") near the mid-point between the fit points. Sometimes this error can be reduced by adding information from the next points out on each end of the interval (can fit a cubic equation). However, think about the general shape for a cubic function. A substantial danger exists in introducing kinks during such a fit and actually decreasing accuracy. I can also use quadratic fits when the first derivative must be continuous at the end points of intervals. If I also need continuity of the second derivative, I fit a series of cubics to the points. During this process, some extra constraints can be imposed to remove or minimize cubic wiggles. This is called CUBIC SPLINE.

Although more computational expense is required for cubic splines, less table points are usually required for the same accuracy as a linear interpolation. You will at times be faced with the need to balance "Lower Cost" linear interpolation with more storage against "Higher Cost" cubic splines needing smaller data tables. It is important to do a careful study for your particular application of the accuracy of both and relative computer time requirements. Many times the actual relative speeds of the implemented methods will surprise you. If you can cut your table size enough with a cubic spline to fit in a computer's "cache memory" (faster than normal memory), you may see a significant increase in speed.

For the following examples, I'm going to deal with interpolation between points in the x,y plane. Given a value of x, the interpolation subroutine will return a corresponding y value based upon interpolation between tabulated pairs of x,y values. The first example is linint1.f demonstrating interpolation when the x values in the table are evenly spaced. Even spacing makes it very easy to locate the bounding points for interpolation, but restricts your ability to adjust point location to optimize the mix of accuracy and table length.

In linint2.f the restriction of even spacing between x points is lifted, and the x and y values of points are stored in the arrays xtab and ytab. The disadvantage to this approach is the need for a search to locate the two tabulated points bounding a given input value of x.

In a large number of applications, the value of x input to the interpolation subroutine changes relatively little from one call to the subroutine to the next. In such cases substantial time can be saved by building a memory of the last table points used into the subroutine (linint3.f). When x values change rapidly and randomly from on call to the next, the best tactic is to use a halving search of the x values in the table to locate bounding points.

Fortran does not guarantee that a variable internal to a subroutine (not in the argument list, a COMMON block, or MODULE) will always retain its value after the subroutine has been exited. To force the value to be available the next time that the subroutine is called you must use a SAVE statement. In linint3.f the location of the lower bounding point for the interpolation is stored in the variable "ilast". The non-executable statement "save ilast" has been included in the subroutine, to insure the value is available the next time the subroutine is called.

The previous examples were special cases where the interpolation table has been built into the subroutine. In applications using a wide variety of interpolated results, it is a good idea to create single interpolation subroutine with the x and y point values and number of points provided through the argument list ( linint4.f ).

I won't go into details of anything beyond simple linear interpolation here. You should simply be aware of the properties and potential pitfalls of a cubic spline. If you are interested in the details of this procedure, look at section 18.6 (p. 485) in your text. When you have a serious application requiring interpolation search Netlib or packages available from the National Center for Atmospheric Research (NCAR) for something appropriate.


Back to the Table of Contents / Home


Written and Maintained by John Mahaffy : jhm@psu.edu