Cubic polynomial is often used as main building block in some or other algorithms of piecewise function approximation. Often we then want to draw the result of the approximation on screen or printer. But we lack primitives to do it.

To draw the cubic polynomial we have to approximate it with the polyline. The maximum distance from polynomial to polyline should not exceed some "tolerance" value.

Let we have a cubic polynomial defined at [x1, x2] interval.

To approximate it with polyline we should do the following:

- Get the deviation polynomial, i.e. the difference between the initial cubic polynomial and the straight line passing through its left and right bound points. This polynomial is either identically equal to zero or has one or two extremum(s) at [x1, x2].
- Evaluate the values of deviation polynomial at extremum points. It its absolute values are lower than the tolerance then the initial cubic polynomial can be approximated with a straight line passing through points (x1, y1) and (x2, y2). Otherwise
- Split the initial interval [x1, x2] on two or three subintervals (depending on extremum’s count) and repeat the procedure recursively from (1) for each of subintervals.

# The code

/// <summary>

/// Approximating Cubic Polynomial with Polyline.

/// </summary>

public static class CubicPolynomialPolylineApproximation

{

/// <summary>

/// Gets the approximation of the polynomial with polyline.

/// </summary>

/// <param name="polynomial">The polynomial.</param>

/// <param name="x1">The abscissas start.</param>

/// <param name="x2">The abscissas stop.</param>

/// <param name="tolerance">The tolerance is the maximum distance from the cubic

/// polynomial to the approximating polyline.</param>

/// <returns></returns>

public static Collection<Point> Approximate(Polynomial polynomial

, double x1, double x2, double tolerance)

{

Debug.Assert(x1 <= x2, "x1 <= x2");

Debug.Assert(polynomial.Order == 3, "polynomial.Order == 3");

Collection<Point> points = new Collection<Point>();

// Get difference between given polynomial and the straight line passing its node points.

Polynomial deviation = DeviationPolynomial(polynomial, x1, x2);

Debug.Assert(deviation.Order == 3, "diff.Order == 3");

if (deviation[0] == 0 && deviation[1] == 0 && deviation[2] == 0 && deviation[3] == 0)

{

points.Add(new Point(x1, polynomial.GetValue(x1)));

points.Add(new Point(x2, polynomial.GetValue(x2)));

return points;

}

// Get previous polynomial first derivative

Polynomial firstDerivative = new Polynomial(new double[] { deviation[1], 2 * deviation[2], 3 * deviation[3] });

// Difference polynomial extremum’s.

// Find first derivative roots.

Complex[] complexRoots = firstDerivative.Solve();

// Get real roots in [x1, x2].

List<double> roots = new List<double>();

foreach (Complex complexRoot in complexRoots)

{

if (complexRoot.Imaginary == 0)

{

double r = complexRoot.Real;

if (r > x1 && r < x2)

roots.Add(r);

}

}

Debug.Assert(roots.Count > 0, "roots.Count > 0");

Debug.Assert(roots.Count <= 2, "roots.Count <= 2");

// Check difference polynomial extremal values.

bool approximates = true;

foreach (double x in roots)

{

if (Math.Abs(deviation.GetValue(x)) > tolerance)

{

approximates = false;

break;

}

}

if (approximates)

{// Approximation is good enough.

points.Add(new Point(x1, polynomial.GetValue(x1)));

points.Add(new Point(x2, polynomial.GetValue(x2)));

return points;

}

if (roots.Count == 2)

{

if (roots[0] == roots[1])

roots.RemoveAt(1);

else if (roots[0] > roots[1])

{// Sort the roots

// Swap roots

double x = roots[0];

roots[0] = roots[1];

roots[1] = x;

}

}

// Add the end abscissas.

roots.Add(x2);

// First subinterval.

Collection<Point> pts = Approximate(polynomial, x1, roots[0], tolerance);

// Copy all points.

foreach (Point pt in pts)

{

points.Add(pt);

}

// The remnant of subintervals.

for (int i = 0; i < roots.Count – 1; ++i)

{

pts = Approximate(polynomial, roots[i], roots[i + 1], tolerance);

// Copy all points but the first one.

for (int j = 1; j < pts.Count; ++j)

{

points.Add(pts[j]);

}

}

return points;

}

/// <summary>

/// Gets the difference between given polynomial and the straight line

/// passing through its node points.

/// </summary>

/// <param name="polynomial">The polynomial.</param>

/// <param name="x1">The abscissas start.</param>

/// <param name="x2">The abscissas stop.</param>

/// <returns></returns>

static Polynomial DeviationPolynomial(Polynomial polynomial, double x1, double x2)

{

double y1 = polynomial.GetValue(x1);

double y2 = polynomial.GetValue(x2);

double a = (y2 – y1) / (x2 – x1);

double b = y1 – a * x1;

if (a != 0)

return polynomial.Subtract(new Polynomial(new double[] { b, a }));

else if (b != 0)

return polynomial.Subtract(new Polynomial(new double[] { b }));

else

return polynomial;

}

}

In the code above I’m using the helper class Polynomial encapsulating operations on polynomials including addition, subtraction, dividing, root finding, etc. It’s ported from "Numerical Recipes in C, 2-nd Edition" book with some additions and bug fixes.

You can download the code sample here:

It is Visual Studio 2008 solution targeted to .NET 3.5. It contains WPF Windows Application project designed to demonstrate drawing of some cubic polynomials. You can select one of the curves from Combo Box at the top of the Window, experiment with tolerance and set appropriate XY Scales.