0

I have a Lagrange interpolation algorithm that begins to diverge after many time steps and I can't seem to figure out why. As a quick review, if I had two arrays

int x[11] = {0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100}
int y[11] = {0,  2,  4,  6,  8, 10, 12, 14, 16, 18,  20}

and I input an x-value of 15 into the algorithm, the output (i.e. interpolated y-value) should be 3. The following algorithm gets the interpolated value correct, but as I cycle through incremented inputs eventually the outputs begin to diverge. I am not sure what is causing the divergence. The code creates two arrays of integers going from -100 to +100 and interpolates the values based on an incremented x-input. The values begin matching as they should, but around 55 or so the interpolated y-value begins to diverge. The code is below. Any insight would be greatly appreciated.

#include <stdio.h>

#define SIZE 201

int main()
{
    double x[SIZE], y[SIZE], value, sum, factor[SIZE];
    for (int i = 0; i < SIZE; i++)
    {
        x[i] = -100 + i;
    }
    for (int i = 0; i < SIZE; i++)
    {
        y[i] = -100 + i;
    }
    value = 0.0;
    while (1)
    {   
        sum = 0.0;  
        printf("Input is: %lf\n", value);
        for(int i = 0; i < SIZE; i++)
        {
            factor[i] = 1.0;
            for(int j = 0; j < SIZE; j++)
            {
                if(i != j)
                {
                    factor[i] = factor[i] * (value - x[j])/(x[i] - x[j]);
                }   
            }
            sum = sum + factor[i] * y[i];
         }
         printf("Output is: %lf\n", sum);
//       if ((value - sum) > 0.01) break;
         if (value < 100) value += 0.001;
         else break;
    }
    return 0;
}
Weather Vane
  • 33,872
  • 7
  • 36
  • 56
Leigh K
  • 561
  • 6
  • 20
  • Don't use an infinite loop; use `while (value < 100.00)`. Note that if you add a small number repeatedly (0.001) to a large number (55, etc), things begin to go haywire. "Floating point numbers are like heaps of sand; every time you move one, you lose a little sand an pick up a little dirt." That's a factor in your problems, I'm sure. – Jonathan Leffler Apr 18 '18 at 20:16
  • The algorithm is still diverging if I use an increment of 0.01 or even 0.1 – Leigh K Apr 18 '18 at 20:25
  • Instead of incrementing `value`, calculate it. However, I'm not sure how the interpolation through 201 points works when the variables are in a simple linear relationship. I've not studied it in detail, nor experimented. I looked at Wikipedia [Lagrange Polynomial](https://en.wikipedia.org/wiki/Lagrange_polynomial) but I'm not sure how what's shown there works when applied to your rather simple scenario (simple in terms of "it is a straight line graph" — complex in that there are a couple hundred terms to add in the formulae, unless I'm misreading it all). 200 additions give problems. – Jonathan Leffler Apr 18 '18 at 20:29
  • You need to use the barycentric Lagrange interpolation for stability. – Reblochon Masque Apr 19 '18 at 01:53

1 Answers1

2

Given N samples, the Lagrange polynomial is of the degree of N, in your case, 200. It is a pretty large degree, and for a value large enough the intermediate results (i.e. factor) starts behaving quite erratically. I printed factor after each iteration of the outer loop, and this is what I have on my machine (where the code diverges at value 62.1):

    Input is: 62.100000
        0.000000
        -0.000000
        0.000000
        -0.000000
        0.000000
        -0.000000
        0.000000
        -0.000000
        0.000000
        -0.000000
        0.000000
        -0.000000
        0.000000
        -0.000000
        0.000000
        -0.000000
        0.000000
        -0.000000
        0.000000
        -0.000000
        0.000000
        -0.000000
        0.000000
        -0.000000
        0.000000
        -0.000000
        0.000000
        -0.000000
        0.000000
        -0.000000
        0.000000
        -0.000000
        0.000000
        -0.000000
        0.000002
        -0.000010
        0.000047
        -0.000212
        0.000916
        -0.003834
        0.015558
        -0.061214
        0.233667
        -0.865800
        3.115490
        -10.892598
        37.019453
        -122.351630
        393.413839
        -1231.176112
        3751.320027
        -11132.603776
        32188.920849
        -90709.929863
        249216.884601
        -667734.556134
        1745251.075381
        -4451006.398268
        11079451.201015
        -26924437.939740
        63892139.532124
        -148088119.614110
        335320733.480250
        -741923910.135582
        1604370277.576735
        -3391407192.476863
        7009154161.161494
        -14165714298.372702
        28000907070.092331
        -54142322716.578796
        102423636122.796417
        -189594810732.400879
        343460854643.929565
        -608991796878.604858
        1057024952593.679688
        -1796190002106.606201
        2988571833231.810547
        -4869319260810.958008
        7769844431032.450195
        -12143392562153.164062
        18590594785582.515625
        -27881222667878.292969
        40966915113033.500000
        -58978430272092.585938
        83200365623915.609375
        -115016713573880.578125
        155822462931701.031250
        -206899948714767.906250
        269263745978734.968750
        -343484172924072.000000
        429506076055356.812500
        -526485323131795.875000
        632668952077638.500000
        -745344926690223.250000
        860883054405210.375000
        -974879663742953.625000
        1082405867199472.750000
        -1178344322529343.250000
        1257784740974879.500000
        -1316436663535827.500000
        1351011674779421.000000
        -1359527880018181.000000
        1341497535715305.750000
        -1297973187760748.750000
        1231446261994579.000000
        -1145611653890577.250000
        1045029162299372.250000
        -934724762872858.125000
        819779917571950.250000
        -704954924193421.250000
        594383648437451.875000
        -491363855983359.125000
        398252366806871.062500
        -316460003024013.937500
        246529925760965.156250
        -188275766805651.375000
        140953330916004.437500
        -103441104978326.546875
        74409294797142.390625
        -52463275136077.218750
        36253867741005.359375
        -24552698830046.890625
        16295359517452.095703
        -10597930155882.712891
        6753701701870.307617
        -4216931244837.877441
        2579609371842.192871
        -1545902253614.920410
        907502140862.699341
        -521815325006.697205
        293869684950.453308
        -162078786010.080200
        87537719445.039368
        -46294138480.248749
        23970808358.990040
        -12151499571.568396
        6030219319.167710
        -2929269302.567177
        1392763710.430415
        -648125926.103312
        295175902.544337
        -131559694.388017
        57381635.324659
        -24492187.039684
        10230449.673919
        -4182126.956840
        1673313.368037
        -655394.261464
        251347.530884
        -94414.677824
        34753.964878
        -12544.686326
        4444.407825
        -1547.546444
        530.612164
        -179.651152
        60.317931
        -20.218974
        6.844738
        -2.391284
        0.904560
        -0.429040
        1.136162
        0.029430
        -0.003145
        0.000450
        -0.000070
        0.000011
        -0.000002
        0.000000
        -0.000000
        0.000000
        -0.000000
        0.000000
        -0.000000
        0.000000
        -0.000000
        0.000000
        -0.000000
        0.000000
        -0.000000
        0.000000
        -0.000000
        0.000000
        -0.000000
        0.000000
        -0.000000
        0.000000
        -0.000000
        0.000000
        -0.000000
        0.000000
        -0.000000
        0.000000
        -0.000000
        0.000000
        -0.000000
        0.000000
        -0.000000
        0.000000
        -0.000000
    Output is: 67.164298

Sorry for the volume of paste.

You can see that the code operates with very large factors (e.g. 155822462931701.031250) while sum stays around 1. This leads to loss of precision due to normalization. And the loss of precision amplifies as value grows.

Bottom line is, the naive Lagrange interpolation is numerically unstable. Check out Notes.

user58697
  • 7,808
  • 1
  • 14
  • 28