Интерполяция назад полиномом Ньютона

[za5]

Всем привет! Дана функция f(x) = 1/(1+25*x*x). Интервал -1..1 .Для нахождения конечных разностей массив узловых точек переворачиваю на изнанку, т.е. нулевой элемент = последнему и т.д. Вот он код: c#

//Нахождение конечных разностей
//ВНИМАНИЕ: нулевой элемент = последний (переназначили)
double[] KonechnR(ref double[] x, ref double[] y, int n)
{
   /* Элементы массива функций в узловых точках будем смещать влево к 
    * нулю, тем самым получая каждый раз в нулевой позиции требуемую
    * конечную разность НУЖНОГО порядка.
    * Найдём все конечные разности и выведем их в KonechnRazn[n],
    * где n - число считаемых конечных разностей */

    double[] KonechnRazn = new double[18];
    KonechnRazn[0] = 1;

    for (int k = 1; k <= n; k++)
    {
        for (int z = 0 ; z <= n - k; z++)
        {
            y[z] = y[z + 1] - y[z];
            if (z == 0) KonechnRazn[k] = y[0];
        }
    }        
    return KonechnRazn;
}    
/* Нахождение значения полинома Ньютона при интерп. назад
 * Исползуем полином Ньютона, записанный с учётом равно-
 * удаленности узлов */
double Newton(double X,double[] x,double[] y,int nOut,double[] KonechnRazn)
{
    double h = Math.Abs(x[1] - x[0]);
    double q = (X - x[0]) / h;
    double result = y[0];
    double proizv = 1;
    double factorial = 1;
    for (int z = 1; z <= nOut; z++)
    {
        factorial *= z;
        proizv *= (q + z - 1);
        result += KonechnRazn[z] * proizv / (factorial);
    }
    return result;
}

Всё достаточно ясно, однако конечный ответ сильно расходится с ожидаемым что бы я не выделывал, может ошибка в вычислении конечных разностей?

[ewert]

za5 в сообщении #359958 писал(а):

    for (int k = 1; k <= n; k++)
    {
        for (int z = 0 ; z <= n - k; z++)
        {
            y[z] = y[z + 1] - y[z];
            if (z == 0) KonechnRazn[k] = y[0];
        }
    }        

Как минимум тут. У Вас по замыслу очередная стартовая конечная разность должна размещаться на очередной освободившейся позиции. Но позиции у Вас освобождаются с конца, а адреса размещения перебираются, наоборот, сначала. Получается перекрытие строчки текущих разностей и строчки накопленных результатов. Измените соответствующим образом индекс в последней строке.

И, кстати, совершенно непонятно, зачем переворачивать исходный массив узловых значений. Куда логически проще и надёжнее перевернуть при необходимости массив полученных разностей.

[za5]

выкинул кусок,

    double Newton(double X,double[] x,double[] y,int nOut)
    {
        double h = Math.Abs(x[1] - x[0]);
        double q = (X - x[0]) / h;
        double result = y[0];
        double proizv = 1;
        double factorial = 1;
        for (int z = 1; z <= nOut; z++)
        {
            for (int j = 0; j <= nOut - z; j++)
            {
               y[j] = y[j + 1] - y[j];
            }
            factorial *= z;
            proizv *= (q + z +1);
            result += y[nOut - z] * proizv / (factorial);
        
        }
        return result;
    }

Но теперь результаты как-то чередуются: первые 12 штук
-587,308459554845
0,185951015770217
-563,058515727851
0,258352684764942
-543,32612466001
0,199822686217371
-527,467143093522
0,114170444855774
-514,897729780225
0,0693863723327772
Если смотреть на четные то близко даже с истиной, но глядя на нечетные напрашивается вывод о перекрывании

[ewert]

Вот тут явная ошибка:

    proizv *= (q + z +1);

(должно быть z-1, если я тоже чего не зевнул). В остальном ошибок вроде не вижу. Ну разве что факториал не было смысла накапливать, но это не ошибка, а нерациональность.

А что это за результаты -- не понял.

[Null]

Эта функция известна тем что при её интерполяции многочленом в точках равномерно расположенных на отрезке полученный многочлен сильно отклоняется от самой функции.

[ewert]

Null в сообщении #359981 писал(а):

Эта функция известна тем что при её интерполяции многочленом в точках равномерно расположенных на отрезке полученный многочлен сильно отклоняется от самой функции.

Сильно-то сильно, но не таким же образом, там большие отклонения вблизи концов и осциллируют вовсе не так, как было написано.

[za5]

Безусловно, там z-1, видимо сохранил на промежуточном этапе, но числа становятся даже больше, решил переделать вот так:

double Newton(double X,double[] x,double[] y,int nOut)
    {
        double h = Math.Abs(x[1] - x[0]);
        double q = (X - x[0]) / h;
        double result = y[0];
        double proizv = 1;
        double factorial = 1;
        for (int z = 1; z <= nOut; z++)
        {
            for (int j = nOut; j >= z; j--)
            {
               y[j] = y[j - 1] - y[j];
            }
            factorial *= z;
            proizv *= (q + z -1);
            result += y[z] * proizv / (factorial);
        } 
        return result;
    }

// Участок делится на 100 частей вот результаты: степень полинома = 14
в точке -0.2 почему то выходит 0.85, а в x= -1 выходит 480
по числам получается лучше... но правильно ли(

[ewert]

za5 в сообщении #359985 писал(а):

x= -1 выходит 480по числам получается лучше... но правильно ли(

как это может быть верно, если отклонение в узле всегда равно нулю?...

Есть гипотеза. Примерно -480 в минус единичке как раз получается, если Вы делите отрезок на 16 частей, а для построения многочлена 14-той степени используете только внутренние узлы (т.е. отбрасываете единичку и минус единичку). Что делать вообще-то нехорошо. А между узлами отклонения на краях тоже большие, но не шибко -- где-то порядка семи.

Да, и ещё не понял: что там за 100 частей?... какие ещё 100?...

[za5]

ewert в сообщении #359988 писал(а):

какие ещё 100?...

при выводе графика подставляю точку интервала [-1;1] в качестве X, и вывожу значения полинома в нескольких точках (100) в файл.

в плане заполнения узловых массивов исправлю и посмотрю ещё правомерность выбора индексов, вот, кстати, заполнение:

//Заполняет массивы узлов (Обратная интерп.)
void AddToArrayBackInterpolation(ref double[] X, ref double[] Y, int n, double xleft, double xright)
{
    double dx = (xright - xleft) / (n++);
    for (int i = 0; i < n; i++)
    {
        X[i] = xright - i * dx;
        Y[i] = Fx(X[i]);
    }
}

-- Чт окт 07, 2010 18:56:25 --

в общем, разумно будет сейчас трассировку сделать...

[za5]

Спасибо всем за помощь, но я решил другим путём пойти, и через разделенные разности:

//Интерполяция назад
double Newton(double point, double[] x, double[] y, int nOut)
{
    double proizv;
    double result = y[0];
    double[,] mass = new double[18, 18];
    for (int i = 0; i <= nOut; i++) mass[i, 0] = y[i];
    for (int k = 1; k <= nOut; k++)
       for (int a = 0; a <= nOut - k; a++)
          mass[a, k] = (mass[a + 1, k - 1] - mass[a, k - 1]) / (x[a + k] - x[a]);
    for (int z = 1; z <= nOut; z++)
    {
        proizv = 1;
        for (int j = 0; j <= z - 1; j++) 
           proizv *= point - x[nOut - j]; 
        result += mass[nOut - z, z] * proizv;
    }
    return result;
}

Всё более-менее прилично выглядит. На этот раз массивы не переворачиваю наизнанку.