Consider a meteorological experiment that measures the temperature of a set of buoys located on a rectangular grid at sea. We can emulate such an experiment by indicating the longitude and latitude of the buoys on a grid of 16 × 16 locations, and random temperatures on them between say 36ºF and 46ºF:
In [1]: import numpy as np, matplotlib.pyplot as plt, \ ...: matplotlib.cm as cm; \ ...: from mpl_toolkits.basemap import Basemap In [2]: map1 = Basemap(projection='ortho', lat_0=20, lon_0=-60, \ ...: resolution='l', area_thresh=1000.0); \ ...: map2 = Basemap(projection='merc', lat_0=20, lon_0=-60, \ ...: resolution='l', area_thresh=1000.0, \ ...: llcrnrlat=0, urcrnrlat=45, \ ...: llcrnrlon=-75, urcrnrlon=-15) In [3]: longitudes = np.linspace(-60, -30, 16); \ ...: latitudes = np.linspace(15, 30, 16); \ ...: lons, lats = np.meshgrid(longitudes, latitudes); \ ...: temperatures = 10. * np.random...
We have three different implementation methodologies to deal with interpolation problems:
A procedural mode that computes a set of data points (in the form of ndarray
with the required dimension) representing the actual solution.
In a few special cases, a functional mode that provides us with numpy
functions representing the solutions.
An object-oriented mode that creates classes for interpolation problems. Different classes have different methods, depending on the operations that the particular kinds of interpolants enjoy. The advantage of this mode is that, through these methods, we can request more information from the solutions: not only evaluation or representation, but also relevant operations like searching for roots, computing derivatives and antiderivatives, error checking, and calculating coefficients and knots.
The choice of mode to represent our interpolants is up to us, depending mostly on how much accuracy we require, and the information/operations that we need afterwards...
Numerically, it is relatively simple to state the approximation problem for the least squares norm. This is the topic of this section.
In the context of linear least squares approximation, it is always possible to reduce the problem to solving a system of linear equations, as the following example shows:
Consider the sine function f(x) = sin(x) in the interval from 0 to 1. We choose as approximants the polynomials of second degree: {a0 + a1x + a2x2}. To compute the values [a0, a1, a2] that minimize this problem, we first form a 3 × 3 matrix containing the pairwise dot products (the integral of the product of two functions) of the basic functions {1, x, x2} in the given interval. Because of the nature of this problem, we obtain a Hilbert matrix of order 3:
[ < 1, 1 > < 1, x > < 1, x^2 > ] [ 1 1/2 1/3 ] [ < x, 1 > < x, x > < x, x^2 > ] = [ 1/2 1/3 1/4 ] [ < x^2, 1...
In this chapter, we have explored two basic problems in the field of approximation theory: interpolation and approximation in the sense of least squares. We learned that there are three different modes to approach solutions to these problems in SciPy:
A procedural mode, that offers quick numerical solutions in the form of ndarrays
.
A functional mode that offers numpy
functions as the output.
An object-oriented mode, with great flexibility through different classes and their methods. We use this mode when we require from our solutions extra information (such as information about roots, coefficients, knots, and errors), or related objects (such as the representation of derivatives or antiderivatives).
We explored in detail all the different implementations for the interpolation coded in the scipy.interpolate
module, and learned in particular that those related to splines are wrappers of several routines in the Fortran library FITPACK
.
In the case of linear approximations in the least squares...