Use LAPACK to solve a linear system and find an interpolating polynomial to construct new points between a series of known data points.

Framework

- Accelerate

## Overview

This article demonstrates how you can generate a continuous curve that passes through a small set of points by computing an *interpolating polynomial*. A polynomial is the sum of a series of terms constructed from variables, coefficients, and exponents (for example, *6x³ + 7x²*, where *6* and *7* and the coefficients, and *x* is the variable); and an interpolating polynomial fills in the gaps between the supplied variables and coeefficients.

For any number of data points, there is a unique interpolating polynomial of order (that is, the largest exponent) which is the number of data points minus one. However, for large numbers of data points, this solution can become numerically unstable.

The image below shows five known points, as white dots, and the values generated by evaluating the found interpolating polynomial, as a red line:

The code in this article determines the polynomial coefficients using a Vandermonde matrix based on the x-components of the known points. The coefficients are the solution to *Ax=b*, where *A* is the Vandermonde matrix and *b* is a vector of the y-components of the known points. You'll use LAPACK to solve *Ax=b*. LAPACK is an acronym for Linear Algebra Package and is a standard software library for numerical linear algebra.

### Generate Known Data

Create an array containing five two-element vectors that describe the known data points between which the code will interpolate.

In a real-world app, you will most likely acquire data points from an external source such a meteorological or financial data source. For this example, specify x-components that are evenly distributed between 0 and 1023, and generate random y-components:

### Create a Vandermonde Matrix

Constuct a Vandermonde matrix where the rows are defined by the elements in a source vector that are successively raised to each integer power up to the the source vector's element count, minus one. For example, in the case of a five-element source vector, *x*, the Vandermonde matrix is of the form:

The Vandermonde matrix used in this article derives from the x-components of the points you're interpolating. For example, given the following points:

The resulting Vandermonde matrix contains the following values:

The following code constructs a Vandermonde matrix from the `points`

array:

### Calculate Coefficients

The coefficients for the polynomial are the solution to *Ax=b*, where *A* is the Vandermonde matrix and *b* is the y-components of the known points. For example, using the matrix created in Create a Vandermonde Matrix, the coefficients are the *x* in the following:

Create the function `solve`

to encapsulate the LAPACK routines that solve *Ax=b*. Note that LAPACK overwrites *b* with the solution vector, *x*:

On return, `coefficients`

contains the polynomial coefficients.

### Use LAPACK to Solve Linear System

Use the LAPACK `dgels`

routine to perform the solve. The `dgels`

name derives from **d**ouble-precision, **ge**neral-matrix, **l**east-**s**quares.

By default, LAPACK expects matrices in column-major format. Specify transpose to support the row-major Vandermonde matrix.

Define constants used by LAPACK routines.

Pass

`-1`

to the`LWORK`

parameter of`dgels_`

to calculate the optimal size for the workspace array.Create the workspace array based on the workspace query result.

Perform the solve by passing the workspace array size to the

`LWORK`

parameter of`dgels_`

.

### Evaluate the Polynomial

The vDSP `evaluate`

function evaluates a polynomial. For example, the following code evaluates a simple polynomial that consists of three variables and three coefficients:

On return, `result`

contains `[18`

by performing the following:

Note that the number of elements returned by the polynomial evaluation is the same as the number of elements in the `variables`

array.

To create an interpolation result that contains 1024 elements, use `ramp(in:`

to create the variables:

On return, `polynomial`

contains 1024 elements with the indices corresponding the the x-components, and the values corresponding to the interpolated y-components: