# Optimal method to numerically integrate a spectral distribution data

I am trying to numerically integrate a spectral distribution data. A spectral distribution function gives the Luminosity per Wavelength of an object (like a galaxy) as a function of the Wavelength. In more mathematical terms,

$$dL(\lambda) = \rho(\lambda) d\lambda$$

where $L$ = Luminosity in ($\lambda+d\lambda$, $\lambda$); $\rho$ = spectral distribution function; $\lambda$ = Wavelength.

The problem is I am not presented with the spectral distribution function per se. What I am given is the value of this function at discrete wavelengths. For example, my data would look like:

==================================
|Wavelength|Spectral Distribution|
|   ...    |         ...         |
|   248    |         500         |
|   300    |         510         |
|   305    |         512         |
|   306    |         518         |
|   ...    |         ...         |
==================================


Please do observe that the wavelength's are not varying in a regular fashion. Up until now I've been just using the following primitive formula to find the total luminosity.

$$\Delta L=\frac{(\lambda_2-\lambda_1)}{2}[\rho_2+\rho_1]$$ $$L = \sum \Delta \lambda$$

Is there any other method that I use which gives an improvement over this method I am using?

EDIT: Here is a plot of $\log_{10}(\rho)$ vs $\lambda$:

-
Could you show a plot of your $\rho(\lambda)$ values? Do you have any plausible functional form for what your $\rho$ looks like (a la Planck, for instance)? – J. M. Jul 31 '11 at 12:07
@J. M. I've added a plot of $\log_{10}(\rho)$ vs $\lambda$. Unfortunately I don't know what the functional form of $\rho$ is. – Bernhard Heijstek Jul 31 '11 at 13:04
Quite noisy, apparently... do you think the usual assumption of only $\rho$ is error-contaminated applies, or do we also have to account for errors in $\lambda$? If the machine determines at which $\lambda$ to sample, then I guess we can assume that the $\lambda$ are error-free... – J. M. Jul 31 '11 at 13:22
@J. M. Actually this is not a real experimental data. This data is from a numerical simulation of a cluster of stars (which in turn uses real data with errors) As you said, the $\lambda$'s can be considered error-free. – Bernhard Heijstek Jul 31 '11 at 13:35
Least-squares may or may not be the best way to determine a $\rho(\lambda)$ that you can subsequently integrate. Does the physics/theory really not suggest any plausible form for $\rho$? I was about to suggest splines until I saw your plot... – J. M. Jul 31 '11 at 13:38

Consider trying to approximate $J = \int_a^b f(x)\, dx$ given noisy values of a smooth function $f$ at $N+1$ equally spaced points $x_i = a + i (b-a)/N$, $i=0 \ldots N$. Thus what you are given is $y_i = f(x_i) + \epsilon_i$ where $\epsilon_i$ are independent random variables with mean 0 and standard deviation $\sigma$. A typical quadrature rule would be of the form $A = \sum_i w_i y_i$ where $\sum_i w_i = b-a$. The error $J - A$ consists of two parts, the truncation error $T = J - \sum_i w_i f(x_i)$ and the random error $R = \sum_i w_i \epsilon_i$. For the Trapezoid Rule (all $w_i = (b-a)/N$) the truncation error is $O(1/N^2)$, and other quadrature rules can make the truncation error go to 0 much faster as $N \to \infty$. On the other hand, the random error has standard deviation $\sigma \sqrt{\sum_i w_i^2}$. which is at least $\sigma (b-a)/\sqrt{N}$ (the minimum being attained for the Trapezoid Rule). This is because $$\sum_i w_i^2 = (b-a)^2/N + \sum_i (w_i - (b-a)/N)^2$$ Thus for large $N$ the random error is much larger than the truncation error, and the Trapezoid Rule (which minimizes the standard deviation of the random error) will work better than other rules that have less truncation error but greater random error.