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.

Introducing the Twisted Plane

This idea is motivated by something that comes up often in the world of surveying. Not always do we have simple straight lines, and not always do we have a simple plane. Here is an example of the situation in plan view:

You may know that three points define a plane and that an arbitrary fourth point, may not be on the same plane. We can unambiguously define the elevation of every point on the outside of our rectangle by "stringing a line" (which is surveyor-speak for straight interpolation). We want the elevations we would get if we pulled a string really tight all the way around the outside of the rectangle. (If you care about the units, call them meters, but they are not really the point of this example.)

What about arbitrary points on the interior of the rectangle?  It is possible to use interpolation schemes to determine the elevation of a point on the interior of a plane such that it will coincide with the plane. But in the above situation, we have four points that are not coplanar.  A common, and highly adaptable, solution to this problem is to break our region up into triangles.  There are two ways we could do this.

Case 1

Notice that in case 1 the elevation at the midpoint of the diagonal line would be 100.125 m.

Case 2

On the other hand, the elevation at the middle of the diagonal line drawn in case 2 would be 100.000 m.

The general case of this sort of procedure - breaking a region up into triangles - is called creating a triangular irregular network (TIN) and is a well established method.  The fact that there may be multiple solutions which are not equivalent does not normally concern us.  We just pick one that works because all we care about is that the water will run.  In a complex plan, it is not uncommon to see the choice of lines given on the plan.  As a surveyor you like to see that as it mean the ambiguity of your job is largely removed.

When it comes to disadvantages, the first thing that is noteworthy about the TIN is that triangles and graders (or their operators) do not like each other very much. Although the grader is a very versatile machine, it still likes following a line. Rectangles are ideal. Look up at case 2 and pretend you're driving a grader, starting from the SW corner and proceeding to the SE corner. As you approach the SE corner, you have a problem: you've got two edges that you are trying to put the blade against. This is frankly infeasible. If this is a very large area (5 miles by 8 miles), we don't care, because the tricky spot is so small compared with the total area.  For a smaller area, this might be significant. So what happens? The watchword is "blend", which basically means the grader operator is supposed to figure it out himself and "make it work."

So, can we model this "blending"? (I don't mean just at the diagonal edge, I want to get rid of that.) Can we figure out what the math looks like that describes the result which the grader (operator) would produce if given four corner points on a small, rectangular area and told to "blend" and "make it work"? If we can do this, then we can apply this same method to areas too big for the grader to cover in just one or two passes. Or apply it to the problem of setting grades on concrete islands that are on such a surface. I call it a twisted plane and it allows us to completely discard the diagonal line.

A twisted plane is a surface formed by using a line (segment) and "sweeping" laterally while changing its angle. The line segment is just like the blade of a grader running from one edge of the rectangle to the opposite side. If the width of our rectangle was exactly the width of the grader blade, it would begin with north end of the blade at elevation 99.50 m and the south end at 100.00 m at the west edge of the area. As the grader proceeded to the east the blade would gradually change its angle from having 0.500 m fall from south to north to having 0.250 m fall from the south to north. (The realities of grading are somewhat swept under the carpet for the sake of simplicity.) Here is an attempt to show what the terrain would look like:

GNU plot of "twisted plane" performed using Maxima

We will now proceed to develop the equation which describes this surface.  Let's start with a fresh diagram that uses variables in it:

For completeness, here are the equations for the elevations of the north, south, east, and west line elevations (we take the south-west corner as the origin):

Now, to determine the elevation at any point on the interior, we can use the west and east lines and do straight interpolation across them like so:

This produces the equation:

If we were to use the north and south lines and interpolate across them instead, we would get:

In both cases the resulting equation is