Exercise 5: Integration

This exercise requires a little independent work; figuring out how to write a function in Maple, and how to implement a function in Fortran. The actual time required, once you know how to do it, is quite small. Spend the rest of the time wisely; ie, have fun!

To help you if you get stuck, there are hints, as usual. If you can do without them, the exercise is more useful, but don't be afraid to use ...

... your bonus points!   Credits: 4/4

NOTE: The next lecture covers material that will be part of the first project, so it is a good idea to come to the lecture.

To extract the exercise files for this weeks exercise, do

    cd ~/ComputerPhysics
    cvs update -d

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

The integrand from the Lecture Notes is available as the function prof.pro in your 5_Integration directory. In this section, the point is to first look at the shape of the integrand and consider what happens with the integrand under the substitution of variables discussed in the Lecture Notes, and then explore how the integration error depends on the number of points.

• Start IDLDE as in previous exercises, and try the following. It is possible, but perhaps just disturbing, to copy and paste from these instructions to the IDL command input; why not just type it in, while actually thinking about what it is that you are doing. Or, if you choose the copy/paste method, use the time you gain to think about the input command, rather than just marching blindly on to the next one :-!

  IDL> n=64                                   ; even number of pts
  IDL> v=100.*(-1.+2.*findgen(n)/(n-1))       ; make equidistant v's
  IDL> plot,v,prof(v)                         ; plot profile as a function of v
  IDL> oplot,v,prof(v),psym=1                 ; use arrow keys + CTRL-a, CTRL-e to add symbols
  IDL> u=atan(100.)*(-1.+2.*findgen(n)/(n-1)) ; construct the u scale
  IDL> v=tan(u)                               ; the corresponding v's
  IDL> plot,v,prof(v),psym=-1                 ; negative psym value -> line + symbols
  IDL> plot,u,prof(v)*(1.+v^2),psym=-1        ; transformed integrand as a function of u
  

• Use IDL to plot the error in the trapezoidal formula as a function of the number of points used (cf. notes). The first command starts a recording of the IDL commands into a journal file error.jou.

  IDL> journal,'error.jou'                                  ; Start recording 
  IDL> plot,/xlog,/ylog,[1,3000],[1e-7,10],xstyle=1,/nodata ; empty frame
  IDL> exact=5.467382771                                    ; from Maple
  IDL> n=2                                                  ; start small
  
  IDL> n=n*2 & v=100.*(-1.+2.*findgen(n)/(n-1))             ; linear v
  IDL> error=abs(trapez(v,prof(v))-exact)/exact             ; relative error
  IDL> print,n,error & plots,n,error,psym=1                 ; print and plot symbol
  

• Use the arrow keys to repeat the last 3 lines again, until n is equal to 2048. Then reset n and do the case with variable substitution:

  IDL> n=2                                                  ; reset n
  
  IDL> n=n*2 & u=atan(100.)*(-1.+2.*findgen(n)/(n-1))       ; linear u
  IDL> v=tan(u)                                             ; non-linear v
  IDL> error=abs(trapez(u,prof(v)*(1.+v^2))-exact)/exact    ; relative error
  IDL> print,n,error & plots,n,error,psym=5                 ; print and plot symbol
  

• Use the arrow keys to repeat the last 4 lines again, until the values for the case with variable substitution are filled out (again up to n=2048). The next statement closes the journal file.

  IDL> journal                                             ; stop the recording
  

• How many points are needed for better than 10^-5 precision in the answer for each of the two methods?

  For the psym=1 case 8 16 32 64 128 256 512 1024 2048 Credits: 1/-1 OK

  For the psym=5 case 8 16 32 64 128 256 512 1024 2048 Credits: 1/-1 OK

• You may now use the journal file (error.jou) to repeat the commands for a PostScript and a png copy. To make the plot more useful edit the journal file removing eventual mistakes and change it such that the plot includes appropriate labels on the x- and y-axis and add your userid as main title on the figure (if you don't know how to, then check the help on the plot command) and a horizontal line showing the 10^-5 threshold. Now make a PostScript file with the plot:

      IDL> set_plot,'ps'              ; PostScript
      IDL> @error.jou                 ; repeat commands from journal file
      IDL> device,/close              ; close the idl.ps file
  

  And a PNG file:

      IDL> set_plot,'z'               ; Set the devise to a Z buffer (incl 3D information)
      IDL> @error.jou                 ; repeat commands from journal file
      IDL> pic=tvrd()                 ; reads the graphical to a byte array
      IDL> set_plot,'x'               ; X windows again
      IDL> write_png,'error.png',pic  ; saves the image in PNG format
  

Starting Maple

To start maple give the shell command "xmaple".  If it is the first time you run xmaple you get a choice btw "worksheet" or "document" -- choose "document".

Evaluating integrals with Maple

To evaluate an indefinite integral, click on "expressions" in the menu list on the left, and then on the first integral symbol.  Fill in an expression instead of f, and hit RETURN. Try this with the expression 1/(1+x²) (use the hat (caret) sign to make the exponent), which should give you arctan(x).

To evaluate a definite integral, choose instead the next integral symbol in the "expressions" menu, and enter values for the limits.  Try for example the interval -100 to 100.

To evaluate a definite integral numerically, just enter the limits in the form of floating point numbers (a number that includes a period), or enclose the expression in evalf(..).

The fact that the integral of 1/(1+x²) is arctan(x) is precisely the reason for the choice of transformation in the example we are looking at in this exercise. The actual (test) profile is similar to 1/(1+x^2) (as a real spectral line profile would also be). Hence the transformation of variables u=arctan(x), or x=tan(u), turns the integrand into a nearly constant function.

• To illustrate this further, lets construct the test function from the Lecture Notes (type these lines into the Maple document and hit RETURN after each):

  f:=3/(1+v²)
   
  y:=1-exp(-f)

• Given that this shows how to write a function of a function in Maple, write an expression for φ that implements the function in the Lecture Notes (please check the formulae in the Lecture Notes (5.pdf). If you get stuck ...

  ... take this hint!   Credits: -1/-1

• Now plot the integrand

  plot(φ,v=-100..100)
   

• and perform the numerical integral:

  evalf(Int(φ,v=-100..100))

What is the value of the integral? 5.46738271 5.46382771 5.46382771 5.467382771 Credits: 2/-2 OK

In this section, we go through the substitution steps in Maple. This is NOT just a repetition of the IDL and Fortran parts; Maple is able to perform analytical differentiation as well as integration, and in a real situation you might use Maple to try to find a function that is sufficiently "similar" to your complicated integrand, and yet is analytically integrable.

• To perform the variable substitution in Maple, first type and execute a cell with

  v:=tan(u)
   

• With Maple, we can find the transformation "weight function" dv/du with analytic differentiation:

  dvdu:=diff(v,u)

  which becomes 1+tan(u)²
   

• Given these pieces, write the expression for the transformed integrand (use evalf(Int( )). If you get stuck ...

  ... use this hint!   Credits: -1/-1

  Compare this with the formula in the Lecture Notes.
   

• As before, it is easy to plot the integrand;

  Plot(φ (1+v^2),u=atan(-100)..atan(100))

  and to evaluate the integral. Check that you get the same answer as before!

• Repeat the plotting from above, but this time plot the graphics to a file that can be included in a report. There are two ways to do this:
  • Click on the graphics window in Maple. Then choose the Plot menu in the top bar of the Maple window and choose one of the few possible choices.
  • Change the plot devise to one of the many available. To do this check the help on plotsetup and look at the examples at the bottom of the page. Use this to create an image of either of the previous plots in png format (it works even it is not listed as one of the available formats) in a file called maple.png. --- CHANGE the x-range to your personal choice before making the plot!

  As a final step, save your Maple document as integrate.mw. The contents of the file should be the result of going through the exercise steps above.

  Check that you can see the maple.png file in your browser or an image viewer program.

  Please upload the maple.png file Locate and upload your maple.png file: Credits: 4/0 OK

  ---

  Fortran exercises

  [about 60 minutes]

  The directory 5_Integration is set up with a main program Integrate.f90 that calls a slightly modified version of QROMB from the Numerical Recipes library. The modified version has two additional parameters; an input parameter that specifies the requested precision, and an output parameter that returns the number of intervals that were actually used. Note that the actual precision may be better or worse than the requested precision; the internal error estimate that is used in QROMB assumes that the integrand is "well-behaved".

  • In a terminal window, change directory to the ComputerPhysics/5_Integration directory. Open the "Makefile" with an editor and choose the appropiate fortran compiler (ifort, gfortran..) by removing the "#" in fromt of the oppropiate "FC = ..." line -- or make your own line with the name of the compiler installed on your system.
  • Compile the code by typing

    lynx:5_Integration> make
    

  • Run the program by typing

    lynx:5_Integration> Integrate.x
    

  • Look at the output and the program Integrate.f90, and figure out how the QROMB results are produced.
  • Add analogous lines with calls to QSIMP and QTRAP and corresponding print out, to make the program show all three cases for each accuracy.
    • The files are there, and the Makefile is OK, so you just have to duplicate a few lines twice and change the function names. As a result, you should now get three lines of results for every accuracy.
    • If you have problems you may use ddd for debugging as you did in the previous exercise.

  • How many intervals did it actually use for a requested precision of 10^-5, and was the precision target actually met?

    For the QTRAP v case 8 16 32 64 128 256 512 1024 2048 Credits: 1/-1 OK QTRAV v target met? yes no Credits: 1/-1 OK

    For the QSIMP v case 8 16 32 64 128 256 512 1024 2048 Credits: 1/-1 OK QSIMP target met? yes no Credits: 1/-1 OK

    For the QROMB v case 8 16 32 64 128 256 512 1024 2048 Credits: 1/-1 OK QROMB target met? yes no Credits: 1/-1 OK

  • Try the same three methods on the variant with variable substitution
    • Here you must write another function (eg profu(u)) in addition to the prof(v) which is already there. Just copy/paste the prof(v) part to the end of Integrate.f90, and change it (a hint is available).

      profu hint!   Credits: -1/-1

    • You also need to change (or duplicate and change) the calls to qtrap, qsimp and qromb so that they use profu, over the right interval in u (rather than v). If in doubt, you may use this hint.

      profu call hint!   Credits: -1/-1

  • How many intervals did it actually use for a requested precision of 10^-5, and was the precision target actually met?

    For the QTRAP u case 8 16 32 64 128 256 512 1024 2048 Credits: 2/-2 OK QTRAP u target met? yes no Credits: 2/-2 OK

    For the QSIMP u case 8 16 32 64 128 256 512 1024 2048 Credits: 2/-2 OK QSIMP u target met? yes no Credits: 2/-2 OK

    For the QROMB u case 8 16 32 64 128 256 512 1024 2048 Credits: 2/-2 OK QROMB u target met? yes no Credits: 2/-2 OK

  ---

  Home work

  [about 2 hours]

  Already before coming to the exercise hours you should have spent an hour or so reading the Chapter 5 of the lecture notes, and the Chapter in Numerical Recipes about Integration. Use additional time after this exercise to look at how the substitution of variables is implemented in IDL, Maple, and Fortran; i.e., take the last and important step in the chain Preparation -> Exercise -> Contemplation.