(Image: Nick84, Wikimedia Commons.)

Solar-weighted Gaussian Quadrature¶

— Notes by Steven G. Johnson, December 2019. If you find this work useful, please cite: Steven G. Johnson, "Accurate solar-power integration: Solar-weighted Gaussian quadrature," arXiv:1912.06870 (2019).

For many problems involving solar energy, e.g. simulating the efficiency of solar cells, it is necessary to compute integrals (also called quadrature in numerical contexts) over wavelength λ of the form:

$$ \int_a^b S(\lambda) f(\lambda) d\lambda $$

for some desired wavelength interval $(a,b)$, where $f(\lambda)$ is application-specific and $S(\lambda)$ is a solar irradiance spectrum: the intensity of sunlight measured on Earth. For example, to evaluate the efficiency of solar cells, typically $f(\lambda) = \lambda A(\lambda)$ where $A(\lambda)$ is the absorption spectrum of the photovoltaic device.

As is shown by the red curve above, S(λ) is an extremely complicated function tabulated experimentally by the American Society for Testing and Materials (ASTM), displaying many sharp dips and oscillations due to atmosopheric absorption by different molecules.

Problem: Lots of expensive function evaluations for accuracy¶

The challenge is that the application-specific f(λ) function is often extremely expensive to evaluate, e.g. for solar cells it involves solving Maxwell's equations (solving a partial differential equation) for each λ.

Even if f(λ) is a smooth function, after multiplying by the solar weight S(λ) one obtains a very badly behaved integrand, requiring a huge number of function evaluations to obtain good accuracy.

Solution: Specialized Gaussian quadrature scheme for S(λ)¶

We can solve this problem very efficiently because the weight function S(λ) is fixed, allowing us to precompute a quadrature scheme that takes the "nastiness" of S(λ) into account analytically, yielding a formula of the form:

$$ \int_a^b S(\lambda) f(\lambda) d\lambda \approx \sum_{i=1}^N w_i f(x_i) $$

where the $N$ points $x_i$ are known as quadrature points and the coefficients $w_i$ are the corresponding quadrature weights. The goal is to find these points and weights so that the integral is accurately approximated for smooth functions f(λ), where the effect of S(λ) is precomputed in $x_i$ and $w_i$.

A well-known method to accomplish this goal is Gaussian quadrature, in which clever theoretical methods are used to construct an $N$-point quadrature scheme that integrates polynomial functions $f(\lambda)$ exactly up to degree $2N-1$. Moreover, Gaussian quadrature converges extremely quickly with increasing $N$ for smooth functions f(λ). All of this is true regardless of the smoothness (or lack thereof) of S(λ), since S(λ) is taken into account in the construction of the quadrature scheme.

Brief review of Gaussian quadrature¶

The construction of Gaussian quadrature schemes is reviewed extensively elsewhere. In brief, however, we first define an inner product $\langle p_1, p_2 \rangle = \int_a^b S(x) p_1(x) p_2(x) dx$ of real polynomials on $(a,b)$. Then we perform a Gram–Schmidt orthogonalization of polynomials $\{1,x,x^2,\ldots\}$ to find orthogonal polynomials $\{q_0,q_1,\ldots\}$ with respect to this weighted inner product. The Gaussian quadrature points $x_i$, amazingly, are the roots of $q_N(x)$. Moreover, as is reviewed in e.g. chapter 37 of Trefethen and Bau, these roots (and also the quadrature weights) can be obtained by solving a tridiagonal $N\times N$ eigenvalue problem whose entries are obtained by a Lanczos algorithm. Essentially, we just need to integrate $2N$ polynomials against S(λ) — since these integrals only need to be done once, and are of cheap functions (polynomials), we can use expensive brute-force methods to perform them.

All of this process, for an arbitrary weight function, is implemented for us in the Julia QuadGK.jl package. I should also mention one other wrinkle: for numerical stability when constructing Gaussian quadrature schemes computationally for arbitrary weight functions, QuadGK applies the Lanczos process to polynomials in the basis $\{T_0,T_1,T_2,\ldots\}$ of Chebyshev polynomials, after rescaling the interval $(a,b)$ to $(-1,1)$, rather than using "textbook" basis $\{1,x,x^2,\ldots\}$.

Below just need to supply QuadGK the solar-irradiance S(λ), and give it a bit of information on how to do the brute-force integrals on interpolated data by breaking it up into intervals. Then it will construct our solar-weighted Gaussian quadrature scheme to any desired order $N$ and for any wavelength interval!

Julia package requirements¶

The following Julia packages must be installed for this notebook to work: QuadGK.jl (numerical integration), Dierckx.jl (cubic splines), PyPlot.jl (plots), and ZipFile.jl (reading zip-compressed data).

Uncomment and execute the following line to install them. (In IJulia, cells starting with ] are commands sent to the Julia package manager.)

Solar irradiance data¶

To begin, let us download and import the ASTM solar irradiance data. We will use the ZipFile.jl package to unpack the .zip file containing the comma-delimited data, and then the DelimitedFiles standard-library package to read the data.

The data consists of 2002 data points for wavelengths (λ) from 280nm to 4000nm:

Let's plot the data to see the "messy" irradiance spectrum:

Since we are planning on integrating this data, it is nice to see how far out in the decaying "tail" we need to go. To do this, let's compute the fraction of the integral:

$$ \mathrm{fraction}(\lambda) = \frac{\int_{280}^\lambda S(\lambda') d\lambda'}{\int_{280}^{4000} S(\lambda') d\lambda'} $$

as a function of λ, superimposed on a normalized spectrum.

(For this plot, we will estimate the integral using a simple Euler approximation.)

From the above plot, we can see that a significant fraction (≈ 1%) of the integral is in the "tail" from 3000–4000nm, so we need the whole spectrum if we want to compute integrals to many decimal places (unless the application's function f(λ) is rapidly decaying too).

Cubic-spline interpolation¶

To compute integrals, we need more than just a set of data points — we need to interpolate this tabulated data so that we can evaluate S(λ) for any λ in (200nm,4000nm).

Since the underlying physical spectrum is a smooth function, unless the data are extremely noisy it makes sense to use a cubic spline interpolation, which interpolates the data by a sequence of cubic polynomials so that the first and second derivatives are also continuous.

However, there is a wrinkle — for Gaussian quadrature, we need to make sure that $S(\lambda) \ge 0$ everywhere (since $S$ must produce a weighted inner product on polynomials). Our data points are all $> 0$, but it is conceivable that a spline interpolant will dip below zero at intermediate points. To ensure that this is not the case, we will fit our cubic spline to $\sqrt{S}$ at the data points, and then square the result.

We will do this using the Dierckx package in Julia, a wrapper around a popular Fortran library that supports irregularly-spaced data like the ASTM solar data.

The "knots" of a spline are the points where the interpolation switches from one polynomial to the next, i.e. points at which some higher derivative is discontinuous. We will come back to this below, since knowing the points of discontinuity is crucial to computing our weight-function integrals quickly and accurately.

Let's plot the spline along with the underlying data in a smaller wavelength interval 400–500nm so that we can see the wiggles:

The spline fit looks reasonable — notice that the rapid oscillations encompass several data points in a row, so they are presumably real features of the data and not simply noise.

Gaussian quadrature construction via QuadGK¶

The QuadGK package does almost all of the work for us in constructing a Gaussian quadrature rule for a given weight function, its x, w = gauss(...) functions that returns arrays of quadrature points x and weights w.

All we need to do is to give it the weight function, in this case our spline-interpolated data. gauss will then laboriously compute a sequence of integrals of this weight function against polynomials.

Although we don't care too much about the speed of this process, since we only need to do it once for a given quadrature rule, it turns out that there is something simple we can do to speed it up. Numerical integration is best for smooth functions, but our spline interpolant is not smooth: it has discontinuities (in some higher derivative) at the "knot" points (≈ data points). To integrate such a piecewise-smooth function accurately, we should simply perform the integral piecewise. The gauss function allows us to pass a quad argument with a custom integration routine, so below we pass a specialized interpquad routine that performs integrals by breaking them up into pieces at the knots. (Remember: we are doing this work only during the construction of the quadrature rule for S(λ), not during the evaluation of the resulting integrals of f(λ).)

The result, below, is a x, w = gaussquad_interpolant(N, X, W, a, b) function that returns the N-point quadrature rule on the interval (a, b) for a weight function interpolated from data points (X[j],W[j]).

Solar quadrature rules¶

Now that we have the short function above, let's apply it to compute a 15-point Gaussian quadrature rule for the S(λ) weight function. This will exactly integrate polynomial f(λ) up to degree 29, and should be good for any smooth functions as long as they aren't oscillating much faster than our 15 points can capture.

Here we calculate it, time how long the calculation takes with Julia's @time macro, and show the quadrature points and weights:

The calculation took about 22 seconds, which is not too bad since we only have to do this once. But actually the calculation is faster than that — most of these 22 seconds are spent by Julia in JIT compiling the computational routines.

Subsequently running the code, now that it has been compiled, is much quicker:

We can also plot the points and weights. Unsurprisingly, they vaguely resemble the solar spectrum itself, with more weight where the solar irradiance is greater:

We'll also compute a 99-point solar-weighted quadrature rule, suitable for integrating much more oscillatory functions.

Notice that it puts few quadrature points around wavelengths like 2500–3000nm where the solar irradiance is almost zero (due to strong absorption peaks). It still puts a lot of data points in the 3000–4000nm range where the irradiance is small but not tiny (≈ 1% of total), because integrating these ranges accurately is necessary if we want numeric integrals that are correct to many decimal places.

Computing the 99-point rule was more expensive, but still took less than 6 seconds. If we need a 999-point rule (to integrate extremely oscillatory f(λ) functions), we could do that too without too much trouble.

Validation and efficiency¶

To illustrate these quadrature rules and see how they compare to simpler methods, let's integrate the test function: $$ f(\lambda) = \sin(2\pi \lambda / 500) $$ which oscillates with a period of 500nm (eight times in the ASTM data's 280–4000nm range).

(This is not too wiggly, but integrating it against S(λ) accurately is very challenging for simpler integration techniques.)

Let's plot it along with the solar spectrum:

Above, the blue dots show the function f(x) sampled at the quadrature points $x_i$ of our 15-point solar-weighted Gaussian quadrature scheme. This is not very many points! Even without the S(λ) weight, extremely simple methods like the trapezoidal rule would not be able to integrate this f(x) accurately with so few points. Capturing the S(λ) accurately as well would be crazy, except that the effect of S(λ) is already exactly taken into account by the solar-weighted quadrature scheme.

Brute force¶

To begin with, we'll just compute the "exact" integral (to nearly machine precision) by the brute-force Gauss–Kronrod quadrature function quadgk provided by the QuadGK package.

We won't be completely naive, however: instead of passing 280, 4000 as the integration limits, we'll pass λ..., which "splats" all of the data points: this causes quadgk to break the integral into 2001 sub-segments where our interpolant Sinterp is smooth, giving an extremely accurate result at the cost of lots of function evaluations. We'll increment a global counter cnt in our integrand routine to keep track of the number of times quadgk evaluates our integrand.

The resulting integral 81.264730... is accurate to nearly 15 significant digits, at the cost of 30,000 function evaluations. That's fine for such a cheap f(λ), but not if evaluating f(λ) required us to solve Maxwell's equations once for each λ!

Solar-weighted Gaussian quadrature¶

Next, we evaluate the same integral using our 15-point and 99-point solar-weighted Gaussian quadrature rules from above, and compare the results to the brute-force "exact" answer above.

Using only 15 points, we now obtain better than 1% accuracy! Using only 99 points, our result is accurate to 13 digits.

Without knowing about the magic of Gaussian quadrature rules with precomputed weights, it might seem unimaginable that one could compute this crazy S(λ)f(λ) integral so accurately with so few points.

Brute force: Trapezoidal rule¶

If you are integrating against tabulated data like we have here, the first technique that might come to mind is the trapezoidal rule. If we have M data points for S(λ), we simply evaluate our integrand S(λ)f(λ) at all M points and interpolate with straight lines in between (adding up the area of trapezoids).

Here, that scheme is actually fairly accurate, to about 4 significant digits, but requires $M=2002$ function evaluations so it is horribly impractical for expensive integrands.

Brute force: Adaptive quadrature¶

Another naive approach might be to simply use an adaptive quadrature method that evaluates the integrand at more and more points until a desired accuracy is achieved. In fact, the quadgk routine provided by the QuadGK package is one such algorithm, based on a Gauss–Kronrod method.

However, to resolve the extremely rapid oscillations of S(λ), especially since the interpolated data is not smooth (the third derivative is discontinuous), adaptive quadrature again takes a huge number of function evaluations.

Here, we see that it takes 5295 evaluations for 1% accuracy, and 12135 evaluations for 0.003% accuracy. Again, this is not practical for expensive integrands.

Convergence rate¶

As mentioned above, the accuracy of Gaussian quadrature typically improves exponentially ("spectral" convergence = faster than power-law) with the number $N$ of quadrature points if $f(\lambda)$ is a smooth function. Let's see that here by plotting the error vs. $N$, both for the $f(\lambda)$ above and for a second function $f_2(\lambda)$ that oscillates $10\times$ more rapidly.

Notice that when the integrand $f(\lambda)$ is rapidly oscillating (or has other "sharp" features), you still need to sample it finely enough to capture the features of $f(\lambda)$. Only after you use a large enough $N$ [to roughly fit $f(\lambda)$ to a polynomial of degree $N$] does the spectral convergence "kick in." (Still, we get 9 digits with just 140 points.)

For this $f_2$, even a 99-point quadrature rule is not enough: sampling $f_2(\lambda)$ at just 99 points clearly leads to aliasing effects in which the function appears to be oscillating at a much lower frequency:

Integrating from 280–1100nm¶

To evaluate the efficiency of solar cells, typically $f(\lambda) = \lambda A(\lambda)$ where $A(\lambda)$ is the absorption spectrum of the photovoltaic device. However, for silicon solar cells the absorption decays extremely (exponentially) rapidly for longer wavelengths, so there is not much point in integrating beyond 1100nm in many applications.

We might therefore want quadrature rules for integrating from 280–1100nm, and this is easy to do with our routines above. Here are 15-point and 99-point quadrature schemes, for example: