ASTRO 735 - Problems

Set 1

Due: Mon., Sep 18, in class

Numerical integration

Let's develop your numerical chops. You are going to need it for this course and for the rest of your astronomical life. Your programming choice for this course is fortran or C. If you do not know how to program in ANY language let me know IMMEDIATELY. If you have never programmed, or don't know C or fortran, consider getting a book on one of these two languages (C is prefered for the non-initiate -- don't be a dinosaur like me). In general, when you are having problems with coding, take a look at the source code to see if it makes sense; refer to your book/manual. If you're stuck, ask your peers, and as a last resort ask me. I get really cranky debugging other people's code.

Once you get the hang of basic programming, check out Numerical Recipes (The Art of Scientific Programming, by Press, Flannery, Teukolosky, and Vettering [Cambridge Press]) -- affectionately known as Numrec. This book has a vaste number of subroutines for doing simple to sophisticated numberical analysis of data. The source code is located in:

/usr/local/recipes_c-kr - C (KR)
/usr/local/recipes_c-ansi - C (ANSI)
/usr/local/recipes_f - FORTRAN

Specific subroutines:

Chapter 4 of Numrec is on numerical integration. Read it closely enough that you understand what is going on. There are a number of schemes for integration. Feel free to choose whatever scheme you'd like. I'd recommend qromb.f which calls polin.f and trapzd.f.

Write your integrator code (a main program which calls the relevant subroutines from Numrec)
Think ahead in your coding: be modular. You would like to write a program which could simply recompile with a new function subroutine.

Once you get things running, I'd also recommend that you consider changing all of the `reals' into double precision. If you don't, round-off errors will bite you sooner than later.

(a) Test your program on a simple analytic function:

integrate f(x) = c * x^n for c = 5, n = 1, over a finite interval x = 1 to 2.

(b) Do something intersting: calculate the comoving emissivity of galaxies

The comoving emissivity, commonly refered to as the luminosity density is simply the integral of all luminosity from all sources in a comoving volume, divided by that volume.

All you have to do is integrate the product of the luminosity function of sources and their luminosity over the appropriate range of luminoisity and you're done. The luminosity function -- Phi(M) or Phi(L) -- will be in units of number per comoving volume (typically Mpc^-3 or Gpc^-3). The sources, in one case of particuarly interest are galaxies. A function which describes the composite galaxy luminosity function reasonably well is called the Schecter function:

Phi(M) = b phi* exp[ -b(alpha+1)(M-M*) - exp(-b(M-M*)) ]

or in luminosity:

Phi(L) dL = phi* (L/L)^alpha exp(-L/L*) d(L/L*)

where b = 0.4 ln(10), phi* is a normalization constant in density, and M* or L* is a normalization constant in absolute magnitude or luminosity, and alpha is a shape parameter which corresponds to the (power-law) slope of the faint end of the luminosity function.

Over-plot Phi(M) for a few values of alpha = -0.75,-1,-1.75 between M*-2.5 < M < M*+10.

Before `going numerical:'

analytically calculate the integral of L Phi(L) dL from L/L* = L'/L* to infinity. Use the result to show that for alpha = -2 the integral does not converge, i.e. the integrated luminosity is infinite.

Note that you should now have a closed-form expression for the integral of L Phi(L) dL. You can check your numerical results here because your closed-form expression can be evaluated via tabulated data (for special cases) or via independent Numrec programs. (This is a hint at what you should get.)

Now it's time to integrate numerically. Integrate L Phi(L) dL from infinite luminosity down to three limiting luminosities equivalent to {M*+2.5, M*+5, M*+10} for alpha = {-0.75, -1, -1.75}. Clearly you can't integrate from infinite luminosity -- what did you do? Is the integral converging? If not, try integrating L Phi(M) dM.

When you're done with this, calculate the three indefinite integral for these values of alpha. Check your result with the close-form expression for the integral in a few cases. Tabulate the results.

last update: Sep 5, 2000

Back to the show.