Use of IF Statements for Interpolation


Assignment :

Examine, compile, and execute examples: inter1.f, inter2.f, and inter3.f If the comments get in the way for you, check inter1-c.f, inter2-c.f and inter3-c.f. Read the Web Page on other Branching statements to prepare for the next lecture.

New Fortran:

WRITE, OPEN, CLOSE

We are now at the stage where we can write some relatively useful programs. However, before we do, I want to introduce a way to improve your flexibility in producing results from a program. Although the PRINT statement always prints results to your screen, you can also write these results to a disk file using the WRITE statement. To save typing, and increase flexibility, WRITE does not refer directly to a file name, but to a unit number. For now we will use a simple form:

write (iunit,*) [variable list]

where the variable list is the same that you would use in the PRINT statement, and "iunit" is an integer constant or integer variable that supplies the unit number. For example to write the values of variables "time" and "distance" to unit number eleven, I would include a statement:

      write(11,*) time, distance
So how do I tell the program what this "11" means in terms of a file name? This is done once any time before the first write to unit 11 with an OPEN statement. If what I really wanted was to put lines containing time-distance pairs, into a file named "range-table", I would have used a statement:

      open(11, file='range-table')
Unless you are doing some fairly fancy file manipulations, it is a good idea to locate this statement at the beginning of your main program, or in a special I/O setup subroutine that is called at the beginning of your main program.

Generally when you open a file, you can let the STOP statement take care of the closing it out neatly. However, in some instances such as my examples' use of another program for graphics, you must tell the computer that you're done with the file, before your program is finished. This is done with the CLOSE statement. When I finish writing to "range-table", I include the statement:

      close(11)
to close down unit 11.

Interpolation and Weighted Transitions

As a demonstration of IF tests, we are going to look at some simple, but useful examples of interpolation. We will return to a slightly different use of interpolation once we know how to use arrays.

In many areas of engineering, information on various physical parameters is obtained from a number of non-overlapping experiments. For example Professor Smith is interested in structure to fluid heat transfer for systems with very low flow velocities. Professor Jones is also interested in the same general topic, but for applications with much higher flow velocities. Both take a lot of data, and both publish curve fits (correlations) to this data. However, Professor Smith only looked at flows with a Reynolds number less than 640, and Professor Jones only looked at data for Reynolds number greater than 2000. I've been given the job of predicting the behavior of a new system where Reynolds numbers can vary from 0 to 3000. How do I use their results?

One solution to this problem was found in the "htcoef" examples, another is to evaluate results in the region of missing data with some form of interpolating function. The simplest interpolation is to take Smith's number for the Nusselt number (measure of heat transfer rate per degree of temperature difference) at Re=640, and Jones' Nusselt number at Re=2000, fit a straight line to the two points, and use that line to get values of Nusselt number for any Reynolds number between 640 and 2000. This is linear interpolation.

One potential disadvantage of linear interpolation is that the resulting curve of Nusselt number versus Reynolds number is not smooth at Re=640 and Re=2000, the derivatives are discontinuous. To impose continuous first derivatives at the two ends of the interpolation region, you need two more parameters to adjust in the fitting function. For a simple polynomial, this takes us from a linear to a cubic fitting function. The idea is to choose the four coefficients in a cubic polynomial (Nu=a3*Re**3+a2*Re**2+a1*Re+a0) so that the polynomial gives the same value and first derivative as Smith's data fit when Re=640, and the same value and first derivative as Jones' data fit when Re=2000. These conditions give four equations to solve for the four unknown coefficients (a0, a1, a2, and a3).

Both linear and cubic interpolations are demonstrated in inter2.f for this problem, and inter1.f for a different problem involving transition between two regions of flat behavior.

A variation on the approach to interpolation is to introduce weight functions that gradually transition from Smith's functional fit to the data to Jones' fit to the data over the transition region of 640<Re<2000. Smith likes a function "f", to give dependence of Nusselt number on Reynolds Number: Nu=f(Re). Jones likes a function "g", to give dependence of Nusselt number on Reynolds Number: Nu=g(Re). I must choose a weight function "w(Re)" that has a value of zero at the low end of the transition region ( Re=640) and increases monotonically to 1 at the high end (Re=2000). Within the transition region, I get the Nusselt number from the equation:

Nu = w(Re)*g(Re) + (1 - w(Re))*f(Re)

The simplest form for "w" is a linear dependence:

w(Re) = (Re - 640.0)/(2000. - 640.)

As with straight end-point interpolation, this introduces discontinuities in the first derivative at the ends of the transition region. Two force continuous derivatives, again we resort to a cubic polynomial. In this case the form of the polynomial can be written without knowledge of "g" or "f". First define an intermediate variable:

      x = (Re - 640.0)/(2000. - 640.)
then the weight function is:

      w = 3*x**2 - 2*x**3
This weighted transition approach is illustrated in inter3.f. I have also made two changes to the form of the programming from inter2.f, to formally reduce the number of arithmetic operations. See if you can find them, but don't get too excited about doing these things. A good optimizing compiler should take care of them for you.

Don't take this exercise too seriously. Smith and Jones are imaginary, and data exists for 640<Re<2000


Back to the Table of Contents / Home


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