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 ...


Subsections  

Getting the exercise files  1 min 
Using IDL to explore the problem  
Integrals with Maple  
Substitution of variables in Maple  
Fortran exercises  
Home Work 
Getting the exercise files 
[about 1 minute] 
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.
Using IDL to explore the problem 
[about 45 minutes] 
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.
IDL> n=64 ; even number of pts IDL> v=100.*(1.+2.*findgen(n)/(n1)) ; 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 + CTRLa, CTRLe to add symbols IDL> u=atan(100.)*(1.+2.*findgen(n)/(n1)) ; 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
IDL> journal,'error.jou' ; Start recording IDL> plot,/xlog,/ylog,[1,3000],[1e7,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)/(n1)) ; 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
IDL> n=2 ; reset n IDL> n=n*2 & u=atan(100.)*(1.+2.*findgen(n)/(n1)) ; linear u IDL> v=tan(u) ; nonlinear 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
IDL> journal ; stop the recording


IDL> set_plot,'ps' ; PostScript IDL> @error.jou ; repeat commands from journal file IDL> device,/close ; close the idl.ps fileAnd 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

Integrals with Maple 
[about 30 minutes] 
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".
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^{2}) (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^{2}) 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.
f:=3/(1+v^{2})
y:=1exp(f)

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

Substitution of variables in Maple 
[about 45 minutes] 
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.
v:=tan(u)
dvdu:=diff(v,u)
which becomes 1+tan(u)^{2 }

Compare this with the formula in the Lecture Notes.
Plot(φ (1+v^2),u=atan(100)..atan(100))
and to evaluate the integral. Check that you get the same answer as before!
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.

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
"wellbehaved".
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.
lynx:5_Integration> make
lynx:5_Integration> Integrate.x
Integrate.f90
,
and figure out how the QROMB results are produced.



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).

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.




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.