Linear Interpolation Made Convenient in Maxima

load(interpol); pts: [[0,5],[1,1.5],[4.4,3.14159],[10,10]]; f(x) := ''(linearinterpol(pts,varname='x)); g(x) := ''(cspline(pts,varname='x)); wxplot2d([f(t)],[t,-5,15]); wxplot2d([g(t)],[t,-5,15]); wxplot2d([f(t)],[t,0,10]); wxplot2d([g(t)],[t,0,10]);

(I usually put the load statement in a separate cell in wxMaxima so that I am not re-executing the load every time I correct my main code and rerun the cell.  Minor point.)

Maxima can turn a set of points into a function for linear interpolation.  For that matter, it can make a cubic spline out of the same points.  I'm more interested in the linear interpolated function today—the other is for your reference.  There are three main points of interest in the above code:

1. linearinterpol() takes points, sorts them by the first ordinate, and produces an expression.
2. linearinterpol() produces expressions which use the characteristic function, rather than taking cases to make it transparent to other Maxima routines.
3. The use of the double single quote operator ('') and the single quote operator (') helps me to define the function f(x) such that x is a formal variable and not embedded as an expression.
4. linearinterpol() produces an expression which extrapolates beyond the specified points (from negative infinity to positive infinity) and it is left to the user's discretion to determine the limits to which this extrapolation is reasonable or applicable.

In a Nutshell

All you need is a list of 2D points.  There is more than one way to specify these, but I find it easiest to use a list of lists.  The main list elements are all lists of 2 representing x and y ordinates.  The function sorts these points according to the x-coordinate before proceeding.  You can specify the variable which will be used in the expression.  Provided that you precede that variable with a single quote operator, you can use a variable that is elsewhere assigned a value.  More on that below.

Transparent Output

The characteristic function, which is often assigned the Greek letter chi ("χ"), returns either 0 or 1.  It is a function of two variables, a set and a value.  If the value is an element of the set, it returns 1, otherwise it returns 0.  In the context of our linear interpolation expression, we have a linear expression for each interval in our set of points and each of these is multiplied by a characteristic function with a different set.  Since the sets used for the characteristic functions form a partition of the real numbers, we won't have more than one of the terms of our overall expression from linearinterpol() evaluate to a non-zero value at a time.  This expression can be integrated, differentiated, or combined with other expressions without having to take special cases.  Take the monkey off your back, let Maxima do it.

Cool, Useful, but Initially Abstruse Operators

The single quote operator (beside the key) is a way to tell Maxima, "Don't evaluate this thing."  If I write x and x is assigned value, then Maxima will do what is normally very helpful to me, namely, make the substitution for me.  If I don't want Maxima to do that, because I actually want the expression x, not the value that x stands for (because it has been assigned a value), I need to write 'x instead. In contrast, the double single quote operator ('') tells Maxima to replace the given expression with the value of the expression. This way, when I use the delayed assignment operator (:=) it does not delay the execution of the linearinterpol() function, only the evaluation of the expression as a function of x is delayed.  You may recall from my previous post that parameters in function declarations are formal.  Even if x is elsewhere assigned a value, f(x) := some_expression_including_x + x, will still work intuitively.

Here is some sample input followed by the output.  This is the best way to clarify in your mind the use of these operators:

x:2;
y:2*x;
z: 2*'x;
d(x) := z;
d(2);
d(3);
h(x) := ''(z);
h(2);
h(3);

(%o46) 2
(%o47) 4
(%o48) 2*x
(%o49) d(x):=z
(%o50) 2*x
(%o51) 2*x
(%o52) h(x):=2*x
(%o53) 4
(%o54) 6

Extrapolation

The expression that linearinterpol() produces does not protect you from yourself.  If will produce a result from negative to positive infinity without telling you that you have gone beyond the boundaries of your initial data set.  Actually, this is as it should be, in my opinion.  The program gives you the freedom to ask all the what-ifs your little heart desires, but it is up to you to discern what is really true to life and applicable to the particular problem you are trying to solve—whether the extrapolation is truly justified by the data is rightly left as a matter for user judgement.

Field Measurement of Quadrilaterals Using Only a Tape Measure

Trapezoids

A common nonrectangular area that needs to me measured in the field is a trapezoid. Sometimes shapes are approximated by a trapezoid, if a rectangle is not considered a reasonable approximation or not easily visualized. The well-known formula for the area of a trapezoid is

                        

where b₁ and b₂ are the lengths of the two parallel sides and h is the distance between them (measured perpendicularly to the parallel sides, of course). It is noteworthy that the (b₁+b₂)/2 can be interpreted as the average of the lengths of the parallel sides. This is the length of a line referred to as the median of the trapezoid. In the field, it is often more convenient/faster to measure the median and the height (h, also called the altitude) than to measure both bases and the height. If the location of the median can be “eye-balled” with sufficient precision for the purposes being met, this may increase productivity by decreasing the number of measurements needed. For more information on the trapezoid and the median, see a description here.

General Case Quadrilateral

There are (at least) two approaches to measuring the area quadrilaterals when nothing is known about the internal angles or “parallelness” of nonadjacent lines. The following diagram displays both methods.

We can determine the area of this shape by making 5 measurements and produce an exact picture, or by taking 3 measurements and getting an approximate area. Both are “area by triangles” methods.

Area by Triangles – Heron Method

If we measure the sides AB, BC, CD, DA, and one of the diagonals (either AC or BD), we can use Heron’s theorem to determine the area of the triangles on either side of the chosen diagonal. We can use the law of cosines and law of sines to determine the angles. Alternatively, a basic AutoCAD drawing using temporary construction-line circles would allow you to determine these results without manual calculations. This would give us a full description of the area. This method also generalizes well to an arbitrary number of sides (but requiring a lot of measurements--don't use this method).

Area by Triangles – By Altitude

On the other hand, we could measure the length of a diagonal (say, AC), leave the measuring tape in place (or replace with a string line) and use a second tape measure to measure from the other vertices (B and D) to the diagonal. As long as the accuracy/precision of the result desired is not too high, we can “eye-ball” perpendicular to the diagonal to get the heights of the triangles to an acceptable accuracy. The area calculation is obvious from there. This method involves fewer measurements, minimal extra equipment, and is quicker to calculate, which may make it more efficient for a quadrilateral. It does not generalize to an n-sided figure as easily as the previous method, but could be done with additional string lines used simultaneously.

Volume Under a Twisted Plane

If you want to know what a twisted plane is, check out my previous post on it.  Basically, a twisted plane is what you get if you take a plane and twist it and the application I have in view is surface modelling.  In that previous post, I developed the equation of this surface as defined by the (x, y, z) coordinates at the 4 corners of a rectangular area (rectangular in plan view).  For review, that equation is

In construction estimation, it is often necessary to measure volumes of material to be excavated or filled in.  To do so, we need to be able to turn elevation and position data into volumes.  A common approach is to find average depths and multiply by the relevant area.  If the depths are different in different areas, the average will need some kind of weighting or the areas treated separately and summed afterwards for a total.  If using a TIN to model a surface, volume calculations are straightforward and fit this idea of using an average depth and multiplying by the area.  (I have the derivation here.)

So, an interesting question that remains, is whether the twisted plane is as well behaved.  The answer is yes, the derivation is straightforward (but tedious to write and is omitted), and the result is

(The result is obtained by integrating z(x, y) for x over [0, length] and then y over [0, width].)