Numerical Integration of Irregular Temperature Readings

Question

I run a website that generates degree days, a specialist form of weather data used for calculations relating to building energy consumption. Without going into unnecessary detail, temperature is a function of time, and degree days are essentially the integral of that function.

At the moment our site calculates the data using an approximation method based on daily average, maximum, and minimum temperatures. But I'm working on improving that method by using finer-grained temperature measurements.

These finer-grained temperature measurements are taken throughout each day, and the recording interval can be anything from 1 minute to several hours - it depends on the weather station making the recordings. For most weather stations, the readings are fairly regular, but they're not completely regular. A weather station might typically record the temperature every half hour, but there will often be extra readings or missing readings, or readings taken at less regular intervals for certain periods.

Initially I've been using the trapezoidal method to numerically integrate the function of temperature against time. It's working pretty well, but I'm wondering if I might be able to improve it.

I'm not a mathematician, and my understanding of numerical integration is only very basic. I understand that Simpson's 1/3 rule and Simpson's 3/8 rule typically work better than the trapezoidal rule when numerically integrating mathematical functions. But real-world temperature readings don't follow an exact mathematical function. Also I understand that Simpson's rules require equal intervals, which my temperature readings don't consistently have.

I wonder if it might make sense to use Simpson's rules to integrate stretches of temperature readings that have 2 or more consecutive time intervals of equal length, and use the trapezoidal rule for stretches of irregular readings. But then I see here (a paper that I don't pretend to understand properly) that the trapezoidal rule can often work better than Simpson's rule for various classes of "rougher" functions. I would guess that outside air-temperature variation would be classed as pretty rough - the temperature jumps up and down throughout the day for all sorts of reasons.

I could probably come up with some way to estimate the effectiveness of various methods, but it's tricky because there's no "right answer" to compare figures against. So I'm trying to figure out what method would make most sense from a theoretical standpoint.

Do you think the trapezoidal rule is likely to be the best approach for me? Or are there other approaches that might make more sense?

Answer 1

Since your temperature readings aren't equispaced, you can't directly apply Simpson's rule; the approach equivalent to this is to construct the parabola that interpolates three consecutive points (i.e., across two panels), and integrate that. The problem with this approach, of course, is that you need to have an odd number of data points (even number of panels) to do this.

You can use the trapezoidal rule, of course, across each panel, but a probably better idea might be to construct a cubic Hermite interpolant for each panel, and then integrate that. An obvious problem is that four conditions are needed to uniquely determine a cubic for each panel (two points and two derivative values), but estimates of derivative values can be constructed from the data such that the piecewise interpolant is locally monotonic; briefly, an piecewise interpolant is locally monotonic if there are no spurious inflection points or extrema within a panel. One approach to estimating derivative values for monotonic interpolation, due to Fritsch and Carlson, is implemented in FORTRAN as the pchip package, and in MATLAB as the function pchip. A more modern approach, and one which may be better is some situations, is due to Steffen. You may have to experiment which of the trapezoidal rule, Fritsch-Carlson, or Steffen would be best for integrating your data.

---

Let me detail the way one uses cubic Hermite interpolation for integrating data:

Each panel is bounded by two points, $(x_i,y_i)$ and $(x_{i+1},y_{i+1})$. The trapezoidal rule consists of constructing the line joining these two points (linear interpolation):

$$f_i(x)=\frac{x_{i+1}-x}{x_{i+1}-x_i}y_i+\frac{x-x_i}{x_{i+1}-x_i}y_{i+1}$$

and integrating that:

$$\int_{x_i}^{x_{i+1}}f_i(x)\mathrm dx=\frac{x_{i+1}-x_i}{2}(y_i+y_{i+1})$$

Integrating with cubic Hermite interpolation can be considered as a further "improvement" of the trapezoidal rule; briefly, in addition to points, one has derivative values $y_i^{\prime}$ and $y_{i+1}^{\prime}$, from which one constructs the cubic

$$g_i(x)=y_i+y_i^{\prime}(x-x_i)+c_i(x-x_i)^2+d_i(x-x_i)^3$$

where

$\displaystyle c_i=\frac{3\frac{y_{i+1}-y_i}{x_{i+1}-x_i}-2y_i^{\prime}-y_{i+1}^{\prime}}{x_{i+1}-x_i}$ and $\displaystyle d_i=\frac{y_i^{\prime}+y_{i+1}^{\prime}-2\frac{y_{i+1}-y_i}{x_{i+1}-x_i}}{(x_{i+1}-x_i)^2}$.

Integrating $g_i(x)$ might look slightly complicated, however, there is a nice expression for the integral:

$$\int_{x_i}^{x_{i+1}}g_i(x)\mathrm dx=\frac{x_{i+1}-x_i}{6}\left(y_i+4g_i\left(\frac{x_{i+1}+x_i}{2}\right)+y_{i+1}\right)$$

whose verification I'll leave up to you.

I'll just add the note that in the case of equispaced data, integration with a piecewise cubic Hermite interpolant is equivalent to integration with the trapezoidal rule plus corrections at the beginning and end.