Exercise 4: Interpolation

In this Exercise you first get a chance to play with various types of one-dimensional (1-D) and two-dimensional (2-D) splines, to get a feeling for how the splines react to motions of the knots. The interactive examples are implemented with a simple form of IDL mouse tracking, and you are invited to inspect and modify the IDL code; there are many circumstances when a quick-and-dirty IDL experiment or animation can be very informative.

The Problem this week consists of a (bug-free this time) set of Fortran-90 procedures capable of iterative table improvement; i.e., seeking out the minimum resolution needed for a 2-D table lookup to given precision. All that is needed is for you to add a 2-D spline routine, patterned after the corresponding IDL routine.

First be sure to pick up your bonus points!   Credits: 5/5

To update the exercise files (If we have made any changes since your first CVS download), do

    cd ~/ComputerPhysics
    cvs update -d

In case of problems, see the CVS update help page.

To start the splines.pro procedure, which demonstrates 1-D splines, change directory to ComputerPhysics/4_Interpolation, and start IDLDE. To allow IDLDE to find all the needed routines that are located in the exercise directory you have to set one and possible both of the following system settings. In the new idlde version 7.0.6 you only need to set Change directory on file open under the IDL -> Editor prefeerences. In 7.0.0 this does not work and you instead need to include the present work directory in the idl search path. This is done by on the IDL -> Paths preference page where you inset the present exercise directory and tick it active before clicking on Apply and OK.

Then type

        IDL> splines

in the command window panel.

The "splines" procedure opens an interactive plot window, with a plot of a function that is initially constant. The function is a spline that passes through a small number of points that may be moved interactively.

• Click near an existing point to move it, or 'drag' points by using click+hold+move.

The type of spline is either one that is exact for 3rd order functions, or a "natural spline" (with vanishing 2nd derivative at the end points), or a local spline, with derivatives computed from the nearest neighbours, or a monotonized version of the local spline.

• Change interpolation method by clicking on the labels at the top of the plot.
  • The difference between natural splines and 3rd order splines is most obvious near the end points. Choose a large y-value at the 3rd point (x=3) , and click back and forth between a natural and a 3rd order spline.
  • Note the similarities and differences between local and monotonic splines; move the point at x=6 slowly up and down, and look at how the slopes at the two neighbouring points change sign (or do not change sign).

Play around to your hearts content, but then spend a little time playing with the IDL source code as well. First, a trivial change:

• One needs to stop the execution to edit, compile, and set breakpoints. To stop, click on QUIT in the plot window.
• Open the splines.pro file in IDLDE (using "Open .." in the File menu). Look at the splines procedure (final page or so) .
• Modify the IDL source code, so the symbol plotted at each point is a triangle instead of a plus (the PSYM parameter to the IDL plot and oplot procedures controls the symbol type).
  • Find the documentation for PSYM in the on-line (? IDL command) or on-line reference guide (a little tricky, look after system variables or plot -- then look for "PSYM" with your browsers "find").
  • Alternatively, just try a few values until you find the triangle :-).

    NOTE: The psym= is used in two places; first to make the initial plot, and then to update the plot (inside the loop) -- you need to change in both places.

    In any case, what is the PSYM value that gives a triangle? 1 2 3 4 5 6 7 8 9 10 Credits: 1/-1 OK

    
    

Now, let's take a look at a characteristic property of 3rd order, non-local splines; the influence of a single function value on the derivatives at all other points. Since the spline problem is linear (the derivatives are the solutions of a linear equation system), it is enough to play with a single point, and look at the effect on all of the derivatives. Already from the plot one can see, that setting a single function value different from the rest causes a "ringing" in the derivative values; the curve is forced to wiggle up and down to pass through all the points with constant values. The derivatives have signs that alternate, and absolute values that decay away from the perturbed point. How fast does the derivative decay? In order to find out:

• QUIT the previous experiment, if you have not already done so.
• Set a break point just somewhere near the oplot command inside the repeat loop (towards the very end of the file -- for example on the cursor command), where the derivative is available in the variable d.
  • Break points may be set by clicking in the grey area to the left of the line in the source code, and choosing Toggle breakpoint and a blue dot appears in the margin.
• Run the code again. When it reaches the breakpoint, it opens a window telling that is wants to change into the debug layout. This rearranges the window layout and switches on a number of new bottoms in the top bar of the Debug sub-window. These are used to control the way to step through the code. Click the small green arrow -- next to a double-bar ("pause") symbol at the top of the "debug" sub-window -- to Resume the run (or type ".c" at the command prompt).
• You may toggle between the run and the debug setup using the IDL and Debug bottoms in the top menubar of the IDLDE window.
• When the program resumes, click near x=3, y=10 in the plot window. The program again stops at the break point, and you have a chance to check the actual values of the derivative.
  • First, just print the values, using "IDL> print,d".
  • Then print d, normalizing with successive values of d(n); e.g.,
    • IDL> print,d/d(4)
    • IDL> print,d/d(5)
    • ...

  How much does the derivative decay between two points? 0.23 0.27 0.33 Credits: 2/-2 OK

  
   

• The decay rate is directly related to the relative size of the elements in the tridiagonal equation system.
  

  Looking through the Section 4 lecture notes, the relative size is 1:3:1 1:4:1 1:6:1 Credits: 2/-2 OK

  
   

• Actually, the decay rate for 1:N:1 elements is given by one of the roots of a 2nd order polynomial (also called the "characteristic polynomial" -- see the PDF copy of the lecture viewgraphs for a discussion) with coefficients a=1, b=N, c=1.

  So, is the decay rate equal to sqrt(3)-2 2-sqrt(5) sqrt(8)-3 Credits: 2/-2 OK

  
   

  How many intervals does it take for the ringing to drop with a factor of about 10^6 10 20 30 Credits: 2/-2 OK

  
   

  
  

In addition to the splines.pro procedure, there is also a similar one for 2-D splines. Start this by

        IDL> splines2d

and set values by clicking in the left of the two windows. This procedure has more to draw, so click, rather than drag.

Study what happens nearby a "diagonal mountain range":

• Set values non-zero (and roughly constant) along a diagonal (and only there -- avoid also the corners).
• Toggle interpolation type by clicking on "TYPE", and observe the chain of pyramids change into a mountain range. Note the small "lakes" that develop along the mountain range.
• Click on "xsurface" to get an xsurface plot that can be rotated interactively.
  • Note that you must exit the xsurface plot (Done button) to add a new point, so add all points first.
• Rotate the view until you can estimate the depth of the "lakes" (RELATIVE to the height of the "mountains").
  • Choosing (in the final xsurface plot) "Shaded surface", "Tools -> Xloadct
    • If it still does not change color you may need to do

          IDL> device, decompose=0
          

      To find out how set this at startup, look at the Absalon IDLDE startup tutorial.

  How deep are the lakes, roughly? 5% 10% 15% Credits: 1/-1 OK

  
  

Write a Fortran function that implements 2-d interpolation (cf. the lecture notes), either by translating the IDL routine spl2d.pro, or from scratch. This should take less than 30 minutes.

WHAT??!!! 30 min?!! You are asking too much of me!! I can do it in 15 minutes -- I know how to change IDL into F-90!! Sure, I could do it in 30 seconds -- I just need to use the idl2f90 filter! Credits: 1/-1 OK

Look at the call to spl2d in the file Lookup_routines.f90 to see what arguments the function is called with. Note that the argument list is not quite the same in Fortran as in IDL; you really need to look at the Fortran call to know which parameters to use, and what to do with them. For example, the Fortran routine is called with one pair of (x,y) points at a time, not with an entire array of values.

You may either choose to use the idl2f90 filter, with

idl2f90 < spl2d.pro > spl2d.f90

Overall structure of the function:

    Module spl2df
    contains

    function spl2d (....,nx,ny)         ! LOOK AT THE CALL (nx and ny are there!)
    implicit none                       ! forces everything to be defined explicitly
    integer::  ...                      ! integer variables
    real:: ...                          ! real variables, some with dimensions

    ! compute indices to the four nearest corners from x,y
    ... [ Here things depend on if you dimension to (nx,ny) or (0:nx-1,0:ny-1)
    ! put limits on the integer addresses
    ... [ use Fortran-s min() and max() functions ]
    ! compute fractions px, py
    ... [ make sure that these allow px larger than 1 when x is outside ]
    ! compute g and h coeffs in all four corners
    ... [ Use cut & paste to take this from the IDL function ]
    ! compute the interpolated value
    ... [ Cut & paste this too from spl2d.pro ]
    ! return the value
    ... [ A couple of changes are needed to do it the Fortran way ]
    end function spld2d

    end module

Need a hint or four?

How do I create this file in the first place, and what's it supposed to be called?   Credits: -1/-1

What was this thing about Fortran and IDL indices being different?   Credits: -1/-1

There is something that I don't understand about these upper and lower end limitations on the indices --- those funny   <  and  >  IDL expressions..   Credits: -1/-1

So, now I have the answer, but how do I return a function value in Fortran?   Credits: -1/-1

So, based on the careful preparations you did before you came (:-), you pretty soon have something you think is a correct function. Well, make sure there are at least no syntax errors, by compiling the function:

• Compile with gfortran (or whatever Fortran compiler you have access to):

  localhost:home> cd ComputerPhysics/4_Interpolation
  localhost:4_Interpolation> gfortran -c spl2d.f90
  

• If there are syntax errors, fix them, and try again
  • NOTE: You don't have to retype the gfortran command; use the up-arrow key to recall it.

• When the function compiles without errors, continue below..

So, how do you know that the program works (apart from a false feeling of confidence :-)? The hint is free of charge: Testing!

There is a small test program, Test.f90, that tests your interpolation routine on simple polynomials, where the result should be exact. In "real life", you should write such test programs your self; they are typically a smaller effort to code than the program piece they can test and thus, with a small extra investment, you obtain a very valuable result; "true" ¹⁾ confidence that the code works.

¹⁾ Well, there are always catches; all a test really proves is that the code works in the test situation. Whether or not one can be confident that it also works in the application depends on how similar the situation is (i.e., on how good -- and bug free -- the test program is :-).

Anyway, in order to compile the subroutine and the test program, and run the test, here is what you do:

• Compile the two together, without the -c option

  localhost:4_Interpolate> gfortran Test.f90 spl2d.f90 
  

• There really should be no errors, since the function compiled OK by itself, but of course, if you misspelled the function name, or something similar...
• The compilation produces an executable file with the default name a.out, so to run it, you do

  localhost:4_Interpolate> ./a.out
  

The program tries successively 1st, 2nd, and 3rd order functions, prints an error estimate followed by OK or NOK. The test covers both the case where all points lie inside (interpolations proper), and the case where some points lie outside (extrapolations). The latter case is there to check if the bounds on the indices are correct.

• If the test comes out with OK for interpolations, but NOK for extrapolations, you probably have a problem with the bounds on the addresses. Free hint:
  • Depending on how you choose to write things, the limits should be from 0 to n-2, or from 1 to n-1.

If it still does not work, you need to use the DDD debugger, but even if it does, you should try it out in this simple situation:

• Recompile the program, but now with the gfortran compiler and a -g option:

  localhost:4_Interpolate> gfortran -g spl2d.f90 Test.f90
  

• If you are using your own laptop, and haven't already done so; install DDD (see the Absalon tutorial page)
• Start DDD from the command line:

   localhost:4_Interpolate> ddd a.out &           <- free hint here!!
  

  • With Windows/CygWin you need to first start the CygWin X server -- use the command startxwin, or find the "XWin" entry in the start menu.

• Click on the "Run" button and watch the same result as before come out in one of the windows
• Click in the left hand side margin of the editor window to make the program stop in your spl2d function
• Now rerun (click "Run"), and watch it stop iin spl2d.
• Step slowly forward (step button) and look at variables in the "variables" window as they are computed (note that you must have passed a line for it to have been executed -- the break point is "at the beginning of the line").
• If you have a problem with the results, use the information you get to locate the problem, and fix it.
  • Also, try to "proof-read" the code; checking against the IDL procedure (if you did not copy it)

The testing phase should leave you with a spl2d that you may be confident works, also in the full program.

To compile and run the iterative table improvement program Lookup.x:

• Run the "make" command.

  localhost:4_Interpolate> make clean
  localhost:4_Interpolate> make
  

  The first command removes the output from the gfortran compilations, to allow the second "make" to work also at fys.ku.dk, where only the ifort compiler is available.

  You will learn more about the make command later, for now it is enough to know that, using a description of the relation between pieces of a program, stored in a file named Makefile, it may be used to simplify the compilation of programs with more than just one or two subroutines. Here, it has been setup to produce the executable file Lookup.x and then run it once.

• The program prints error estimates for each iteration, where the number of points in each dimension increases from N to (2*N-1) in each iteration.

  How many points are needed for less than 0.0001 relative error in the table? 19 37 73 145 289 Credits: 10/-2 OK

  
  

• After that, the program prompts for temperature and log pressure pairs, and prints the exact and interpolated result. The Makefile gives it one pair of input values.
• To run the program again, do

  localhost:4_Interpolate> ./Lookup.x
  

  • Here you have a chance to try other values, for example outside of the actual table, to see how good or bad the interpolation is when it turns to an extrapolation.
  • To get out of the loop, give Ctrl-d (or Ctrl-c) instead of the input values.
  
  
• If you encounter numerical problems in this phase, repeat the debugging procedure, and figure out what is wrong.
  • Even if there are no errors, you may want to play with the debugger on this program, that has more levels, and hence is more interesting to debug.

After you are done with the exercise, clean up the directory, to save space. Remove *.o, *.x, idl.ps, and a.out files. A convenient way of doing this is built into the Makefile; the target clean does the cleanup, so type

     localhost:4_Interpolate> make clean

Upload the spl2d.f90 file below.

If you don't feel up to date with Fortran, then spend time reading in your Fortran book, about declarations, subroutines and functions, and about arrays and array indices. You should also follow up on the last lecture, by browsing the Chapter about interpolation in Numerical Recipes. Note, however, that the spline discussion there is unnecessarily complicated--consider it cursory material. As a preparation for the next Lecture, it is also a good idea to browse the Chapter in Numerical Recipes about numerical integration (quadrature).