Click here to Skip to main content
15,868,141 members
Articles / Programming Languages / C#

Interpolation from Polynomial to Natural Splines

Rate me:
Please Sign up or sign in to vote.
4.96/5 (48 votes)
25 Dec 2016CPOL12 min read 54.4K   5.7K   78   20
Interpolation from polynomial to natural splines

Introduction

The question about interpolation recently came to me. So I started working on this topic.

I implemented the Polynomial, Lagrange, Newton and natural spline algorithm and started to compare each one to another. And to my great surprise, I found the fancy spline algorithm not to be the best solution in any case. With certain wave shapes, even a polynomial interpolation works better than a spline interpolation.

But first, some words to explain to the algorithms...

Polynomial Interpolation

The polynomial interpolation algorithm builds for n supporting points a polynomial of the degree n like:

Image 1

Where x and y are the coordinates of one supporting point. For n supporting points, we get n such equations for x1, y1 to xn, yn. So the algorithm basically has to set up the equation matrix of n*n and solve this by a Gauss algorithm. That gives the parameters a0 to an. With this parameters for any xp the corresponding yp value can be calculated.

The polynomial interpolation is the easiest algorithm to be implemented of the 4. But it gets to its limits regarding accuracy quite soon. If the deltaX between the supporting points is too small or too big, the Gaussian algorithm gets problems with the constellation of the matrix equation already with 10 supporting points. It can help to scale the deltaX up or down to get a deltaX close to 1 only to solve the matrix equation. I made most of the interpolations with deltaX from 0.1 to 0.4. That worked quite well. But the problem remains. With more than 10 supporting points, the polynomial interpolation is not reliable any more. As can be seen in the following graphic:

Image 2

Trial of a Polynomial interpolation with f=1/x^2 and 20 supporting points. This phenomenon started at 16 supporting points already.

But this is not the end.

In the equation matrix, we start with:

Image 3

go on to the equation of y2, the y3 and so on. That means in the first line, at the very top, we have the smallest values because of x1 and increasing values further down and Gauss does not like that. Filling the equation matrix upside down like...

Image 4 ...

Image 5

Image 6

Image 7

...improves the accuracy of the polynomial algorithm quite a bit. Like this, the 20 supporting points can be interpolated correctly with a deltaX of 0.5.

C#
    // fill matrix upside down
    for (i = 0; i < order; i++)
    {
        m.y[i] = y_k[order - 1 - i];
        for (k = 0; k < order; k++)
        {
            m.a[i, k] = Math.Pow(x_k[order - 1 - i], k);
        }
    }
...

Just filling the matrix like this instead of:

C#
    // fill matrix reular
    for (i = 0; i < order; i++)
    {
        m.y[i] = y_k[i];
        for (k = 0; k < order; k++)
        {
            m.a[i, k] = Math.Pow(x_k[i], k);
        }
    }
...

Improves the accuracy quite a bit.

Lagrange Interpolation

Lagrange builds a polynomial for each supporting value including the xp value for the interpolation. For n supporting values, that makes n polynomials L1..Ln.

Image 8

xp is the interpolation variable to which the value yp shall be calculated.

yp becomes now the sum of all the Lk* yk terms.

Image 9

In the polynomial for Lk can be seen that if xp = xk the polynomial for Lk becomes 1 and yp becomes yk. and the closer xp gets to one of the supporting points, the more weight this point gets.

The interpolation for one interpolation point is made in one loop including the calculation of all the Lk elements and adding the multiplication of Lk * yk together:

C#
    for (i = 0; i < order; i++)
    {
        nominator = 1.0;
        denominator = 1.0;
        /* calculate the Lk elements */
        for (k = 0; k < order; k++)
        {
            if (k != i)
            {
                /* nominator and denominator for one L element */
                nominator = nominator * (xp - x_k[k]);
                denominator = denominator * (x_k[i] - x_k[k]);
            }
        }
        /* put every thing thogether to yp*/
        if (denominator != 0.0)
            yp = yp + y_k[i] * nominator / denominator;
        else
            yp = 0.0;
    }var i = 0;
...

This loop calculates one yp interpolation value.

The Lagrange Algorithm is much more reliable than the Polynomial if the number of supporting points exceeds 10. It also gets its problems with small a deltaX. But 60 supporting points with a deltaX = 0.1 can be interpolated without problems.

Newton Interpolation

Newton uses another polynomial for the interpolation of yp.

Image 10

Image 11

The coefficients b0 to bn he calculates like:

Image 12

Image 13

Image 14

Image 15

Image 16

In this form, the calculation is not too obvious. But if we put the |XkXk-1|… terms into a table (here for n = 4), things become clearer.

Image 17

 

We first have to calculate the |XkXk-1| terms. From this the |XkXk-1Xk-2|…terms and so on. Even though we only need the value on the very top of each column for b0 to b4, we have to calculate all the elements. It seems quite complicated. But in a recursive function, it becomes quite simple.

C#
void CalcElements(double[] x, int order, int step)
{
    int i;
    double[] xx;
    if (order >= 1)
    {
        xx = new double[order];
        for (i = 0; i < order-1; i++)
        {
            xx[i] = (x[i + 1] - x[i]) / (x_k[i + step] - x_k[i]);
        }
        b[step - 1] = x[0];
        CalcElements(xx, order - 1, step + 1);
    }
}

This function is called with the array containing the supporting points, the number of supporting points as the order and it starts with the first step.

Now one yp interpolation value can be calculated like:

C#
for (i = 1; i < order; i++)
{
    tempYp = b[i];
    for (k = 0; k < i; k++)
    {
        tempYp = tempYp * (xp - x_k[k]);
    }
    yp = yp + tempYp;
}

The Newton Algorithm is quite insensitive to a small deltaX and works really fine till 60 supporting points with a deltaX=0.001. I think it is the best of these 3 algorithms.

For more complex problems, the interpolation with natural splines is meant.

Natural Splines

For the spline interpolation, one interpolation function is calculated for each interval between two supporting points.

Image 18

For each of these intervals, one cubic polynomial is calculated like:

Image 19

for one interval starting at xi and ending at xi+1 and x as the interpolation variable.

To make sure that the wave shapes of all these intervals fit together, there are some conditions to be fulfilled:

Image 20 for i = 1 to n-1

Image 21 for i = 2 to n-1

Image 22 for i = 2 to n-1

Image 23 for i = 2 to n-1

From this conditions and the substitution hi = xi+1 - xi the parameters a, b, c and d can be derived:

From condition a, we get ai = yi.

From b, we get:

Image 24

And c yields:

Image 25

These two equations put together lead to a matrix equation of the order n-2 for ci with the matrix...

 

Image 26

...and the solution vector...

Image 27

This matrix equation solved by a Gauss algorithm yields c2 to cn-1. The first and the last c both are set as 0.

Now there is only di left and this we get from condition d:

Image 28

With these parameters, each of the intervals between the supporting points can be interpolated now and at the end, the pieces have to be put together to one graph.

The natural spline algorithm works fine till 100 supporting points. Then it gets its problems with accuracy. But it is surely the most complicated algorithm to be implemented.

Speed improvement for Splines

To calculate the h parameters we have to solve a matrix equation of following structure:

 

x1

x2

x3

x4

 

a1

b1

 

 

= d1

c2

a2

b2

 

= d2

 

c3

a3

b3

= d3

 

 

c4

a4

= d4

Matrix equation for n = 4

This is a tridiagonal matrix (at least as long as we don’t want to interpolate a periodic function). This special case of a matrix equation can be solved much quicker than by using Gauss by transforming the A matrix of the equation Ax = D into a left triangular matrix and a right triangular matrix as A = LR.

a1

b1

 

 

 

1

 

 

 

 

m1

r1

 

 

c2

a2

b2

 

 

l1

1

 

 

 

 

m2

r2

 

 

c3

a3

b3

=

 

l2

1

 

*

 

 

m3

r3

 

 

c4

a4

 

 

 

l3

1

 

 

 

 

m4

 

Comparison of the parameters gives:

 

a1 = m1

b1 = r1

c2 = l1 * m1

a2 = l1 * r1 + m1

b2 = r2

c3 = l2 * m2

a3 = l2 * r2 + m2

b3 = r3

c4 = l3 * m3

a4 = l3 * r3 + m3

b4 = r4

 

And from this:

            li = ci+1 / mi

            mi+1 = ai+1 – li * bi

            for i = 1 to n-1

 

To solve this equation we first extend the equation LRx = D by y and say LRxy = Dy (y is not the solution vector of the very top equation). This can be separated to Ly = D and Rx = y.

Ly = D looks like

1

 

 

 

 

y1

 

d1

l1

1

 

 

 

y2

 

d2

 

l2

1

 

*

y3

=

d3

 

 

l3

1

 

y4

 

d4

 

And resolving Ly = D gives

 

y1 = d1

  • y1 = d1

y1 * l1 + y2 = d2

  • y2 = d2 - y1 * l1

y2 * l2 + y3 = d3

  • y3 = d3 – y2 * l2

y3 * l3 + y4 = d4

  • y4 = d4 – y3 * l3

 

Rx = y looks like

m1

r1

 

 

 

x1

 

y1

 

m2

r2

 

 

x2

 

y2

 

 

m3

r3

*

x3

=

y3

 

 

 

m4

 

x4

 

y4

 

And resolving Rx = y gives

x1 * m1 + x2 * r1 = y1

  • x1 = (y1 – x2 * r1) / m1

x2 * m2 + x3* r2 = y2

  • x2 = (y2 – x3 * r2) / m2

x3 * m3 + x4 * r3 = y3

  • x3 = (y3 – x4 * r3) / m3

x4 * m4 = y4

  • x4 = y4 / m4

 

With this we basically have to implement 3 loops to solve the tridiagonal Matrix equation:

 

First loop:

m1 = a1
For i = 1 to n-1
            li = ci+1 / mi
            mi+1 = ai+1 – li * bi

 

Second loop:

y1 = d1
for i = 2 to n
            yi = di – li-1 * yi-1

 

Third loop:

xn = yn / mn
For I = n-1 to 1
            xi = (yi - bi * xi+1)/mi

 

Now the first and second loop can be implemented in one if we separate the first term of the first loop and the last term of the second loop

 

m1 = a1
l1 = c2 / m1
m2 = a2 – l1 * b1

For i = 2 to n-1
            li = ci+1 / mi
            mi+1 = ai+1 – li * bi
            yi = di – li-1 * yi-1

yn = dn – ln-1 * yn-1

 

So there remain 2 loops to solve the entire matrix equation and that takes a bit more than half as much calculation time as a Gaussian algorithm.

I impelemnted this algorythm in the sample project Splien_LR ind the function SolveLR() :-)

 

Comparison

 

The comparison of these 4 algorithms was a little surprising. The Polynomial, the Lagrange and Newton algorithm tend to overshoot in certain conditions. That’s written everywhere. And that’s true. As long as the supporting point values jump up and down, the interpolated wave does it even more. As can be seen at a sample shape of 6 supporting points:

y

0.5

11.7

14.8

4.0

2.2

0.2

x

0.1

0.2

0.3

0.4

0.5

0.6

Image 29

All 3 algorithms show some overshooting.

The spline interpolation yields a much better result:

Image 30

A shape I found in a book was 4*exp-(x-0.4)^2.

y

0.0005

0.073

1.472

4.0

1.472

0.073

0.0005

x

0.1

0.2

0.3

0.4

0.5

0.6

0.7

Seems really to be too much for the Polynomial, Newton and Lagrange algorithm:

Image 31

And the spline algorithm handles it much better:

Image 32

The spline curve tends to approach the supporting points more straight and to bend in a smaller radius. That looks quite good. But it is not the best solution for all cases.

A wave shape like 1/x with the points:

y

1.0

0.5

0.3333

0.25

0.2

0.1666

0.143

0.125

x

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

Becomes quite different:

Image 33

Here the Polynomial, Newton and Lagrange algorithm make a smooth, round shape. And what does the spline do?

Image 34

Great surprise: That’s more like a banana shape.

So let’s go a step further and look at 1/x^2.

y

1.0

0.25

0.111

0.0625

0.04

0.0277

0.0204

0.0156

x

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

Image 35

Again a smooth shape from Polynomial, Newton and Lagrange.

And the spline curve:

Image 36

O.k. it’s not like a banana anymore. :)

Using the Code

I implemented all 4 algorithms as much the same as possible. There is a global const order” defining the number of supporting points. This number can be set to an arbitrary value. There are 10 input fields for the y value and 10 for the x value of the supporting points. As long as the order is smaller or equal to 10, these fields are editable and the supporting points can be manipulated. If orders bigger than 10 should be tested, the supporting values have to be set in the software (I kept the code as simple as possible). The supporting points are implemented in the arrays x_k and y_k and the interpolated wave is put in the arrays xp and yp containing “datapoints” elements. The global const “deltaX“ might be clear. But it's not a must to use it. The x_k values can be defined free.

The wave shape is drawn in the panel pGraph by the function DrawGraph(int maxPoints, double minX, double maxX, bool doClear). This function scales the graphic to fit into the panel and clears the panel before drawing if doClear is true.

The calculation is done in the button1_Click event.

Points of Interest

It really depends quite a lot on the wave shape that should be interpolation which algorithm is best and it might be wise to compare the algorithms and to have a look at the interpolated shape. For small numbers of supporting points, the simple Polynomial algorithm might be the best solution. But already at the amount of 10 supporting points, a Lagrange or Newton Algorithm should be used whereas the Newton should be preferred. The spline algorithm is quite fancy but it’s not always worthwhile to use it as it is quite complicated to be implemented and not too clear to be debugged (just in case :-)).

 

For more detailed descriptions visit www.mosismath.com :-)

 

License

This article, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)


Written By
Tester / Quality Assurance Annax Switzerland AG
Switzerland Switzerland
Computers are very straight... They always do exactly what we tell them to do... Only, much too often what we tell them to do is not really what we want them to do Smile | :)

Writing Software is one of the most creative tings one can do. I have been doing this for more than ten years now and still having a lot of fun with it. Besides doing software for HMI's on C# for business, I enjoy very much to implement interesting algorithms and analyse the mathematics they are based on in my leisure time Smile | :)

For more detailed descriptions and math visit me on my own page

www.mosismath.com

Comments and Discussions

 
QuestionVery nice! Pin
Sing Abend5-Jun-18 8:19
professionalSing Abend5-Jun-18 8:19 
GeneralInterPolation usage in FEM analysis Pin
Member 43208441-Jan-17 8:13
Member 43208441-Jan-17 8:13 
PraiseGreat Pin
George196125-Dec-16 10:55
George196125-Dec-16 10:55 
QuestionCool Pin
puromtec116-Dec-16 1:34
puromtec116-Dec-16 1:34 
QuestionSpline interpolation of more than 2 dimensional data arrays Pin
Member 1222807128-Dec-15 22:57
Member 1222807128-Dec-15 22:57 
QuestionDont you think splines are dead? Pin
Alexey KK28-Dec-15 11:03
professionalAlexey KK28-Dec-15 11:03 
Questionspline interpolation Pin
knott9928-Dec-15 9:21
knott9928-Dec-15 9:21 
PraiseMy Vote of 5 Pin
JoeSoMD28-Dec-15 6:51
JoeSoMD28-Dec-15 6:51 
GeneralRe: My Vote of 5 Pin
Mosi_6228-Dec-15 23:12
professionalMosi_6228-Dec-15 23:12 
QuestionIs Lagrange algorithm more stable? Pin
RAND 45586627-Dec-15 4:14
RAND 45586627-Dec-15 4:14 
AnswerRe: Is Lagrange algorithm more stable? Pin
Mosi_6228-Dec-15 23:05
professionalMosi_6228-Dec-15 23:05 
AnswerThank! Pin
RAND 45586629-Dec-15 5:52
RAND 45586629-Dec-15 5:52 
QuestionThe Normal Curve Pin
dpminusa6-Mar-15 14:36
dpminusa6-Mar-15 14:36 
AnswerRe: The Normal Curve Pin
Mosi_628-Mar-15 1:35
professionalMosi_628-Mar-15 1:35 
GeneralRe: The Normal Curve Pin
dpminusa9-Mar-15 16:30
dpminusa9-Mar-15 16:30 
GeneralRe: The Normal Curve Pin
Mosi_6212-Mar-15 9:20
professionalMosi_6212-Mar-15 9:20 
GeneralRe: The Normal Curve Pin
dpminusa12-Mar-15 18:59
dpminusa12-Mar-15 18:59 
GeneralRe: The Normal Curve Pin
Mosi_6213-Mar-15 9:10
professionalMosi_6213-Mar-15 9:10 
GeneralMy vote of 5 Pin
loyal ginger6-Mar-15 10:46
loyal ginger6-Mar-15 10:46 
GeneralRe: My vote of 5 Pin
Mosi_6215-Mar-15 8:19
professionalMosi_6215-Mar-15 8:19 

General General    News News    Suggestion Suggestion    Question Question    Bug Bug    Answer Answer    Joke Joke    Praise Praise    Rant Rant    Admin Admin   

Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.