c - C中的拉格朗日插值算法在许多步骤后发散

Question

我有一个拉格朗日插值算法，它在许多时间步后开始发散，我似乎无法弄清楚原因。快速回顾一下，如果我有两个数组

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}

并且我在算法中输入 15 的 x 值，输出(即插值 y 值)应该是 3。以下算法得到正确的插值，但是当我循环递增输入时，最终输出开始发散。我不确定是什么导致了分歧。该代码创建了两个从 -100 到 +100 的整数数组，并根据递增的 x 输入对值进行插值。这些值开始匹配，但在大约 55 左右的插值 y 值开始发散。代码如下。任何见解将不胜感激。

#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;
}

Answer 1

给定 N 个样本，拉格朗日多项式的次数为 N，在您的例子中为 200。这是一个相当大的次数，对于 value 足够大，中间结果(即 factor)开始表现得非常不稳定。我在外循环的每次迭代后都打印了 factor，这就是我机器上的内容(代码在 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

抱歉糊状物的体积。

您可以看到代码使用非常大的因子(例如155822462931701.031250)运行，而sum保持在1左右。这导致损失由于归一化而导致的精度。随着 value 的增长，精度的损失会增加。

底线是，朴素的拉格朗日插值在数值上是不稳定的。查看 Notes。