VPF::Методы интерполяции - Форум программистов

Shanks

Дата 7.4.2011, 12:14 (ссылка)

(нет голосов)

Загрузка ...

Новичок

Профиль
Группа: Участник
Сообщений: 4
Регистрация: 7.4.2011

Репутация: нет
Всего: нет

Привет всем. Взял алгоритм на с++ с этого форума интерполяции Лагранжа.
Мне необходимо по точкам построить график ну и знать примерное значение в любой точке диапазона. Подозреваю, что из за логарифмического масштаба по оси х алгоритм глючит. Как можно обойти данную проблему ? Возможно ли использовать другие методы ?
Excel ведь как-то строит нормальный график :(
Данные, что я использую есть в коде

Код


#include <cstdlib>
#include <iostream>
 
using namespace std;
 
 
//////////////////////////////////////////////////////////////////////////////
//
//  Interpolation (using Lagrange polynom)
//  (c) Johna Smith, 1996
//
//  Method description:
//    Given: several points of the unknown function
//    We should find values of this function between given points and
//    approximate  function by polynom.
//    Interpolation polynom calculated by the following formula:
//             N       (x-x0)...(x-x(i-1))(x-x(i+1))...(x-xn)
//       L(x)=SUM y(i)----------------------------------------
//            i=0     (xi-x0)..(xi-x(i-1))(xi-x(i+1))..(xi-xn)
//    To calculate values of function between given points we should
//    calculate value of interpolation polynom in these points.
//
//////////////////////////////////////////////////////////////////////////////
 
#include <stdio.h>
 
#define N  8
 
float x[N]={1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 3e9};
float y[N]={-17, -21, -43, -75, -96, -120, -142, -174};
float xvalues[N-1],f[N-1];
int i,j,z;
float l,k;
float pow[N],temp[N],tempaux[N];
 
int main(int argc, char *argv[])
{
  // initializing variables
  for(z=0;z<N-1;z++)
  {
    xvalues[z]=(x[z+1]+x[z])/2; // calculating function between given points
    f[z]=0;
  }
  for(z=0;z<N;z++) pow[z]=0;
 
  // calculating function values
  for(z=0;z<N-1;z++)
    for(i=0;i<N;i++)
    {
      l=1;
      for(j=0;j<N;j++)
        if(i!=j) l*=(xvalues[z]-x[j])/(x[i]-x[j]);
      l*=y[i];
      f[z]+=l;
    }
 
  // Determining interpolation polynom
  for(z=0;z<N;z++)
  {
    k=y[z];
    for(i=0;i<N;i++)
      if(i!=z) k/=(x[z]-x[i]);
    for(i=0;i<N;i++) temp[i]=0;
    temp[0]=k;
    for(i=0;i<N;i++)
      if(i!=z)
      {
        for(j=0;j<N-1;j++) tempaux[j+1]=temp[j];
        tempaux[0]=0;
        for(j=0;j<N;j++) tempaux[j]-=temp[j]*x[i];
        for(j=0;j<N;j++) temp[j]=tempaux[j];
      }
    for(i=0;i<N;i++) pow[i]+=temp[i];
  }
 
  // printing results
  for(z=0;z<N-1;z++) printf("F(%f)=%f\n",xvalues[z],f[z]);
  printf("\nPolynom: F(X)=");
  for(z=N-1;z>0;z--) printf("%fX^%d+\n",pow[z],z);
  printf("%f\n",pow[0]);
  system ("pause");
}

Выводит эти значения и уже на третьем отрезке наблюдается значительное отклонение

Код


F(55.000000)=-19.044754
F(550.000000)=-35.872349
F(5500.000000)=84.443520
F(55000.000000)=-98900.929688
F(550000.000000)=556564096.000000
F(5500000.000000)=-32399998058496.000000
F(1504999936.000000)=10278410099892132000000000000.000000
 
Polynom: F(X)=-0.000000X^7+
0.000000X^6+
-0.000000X^5+
0.000000X^4+
-0.000000X^3+
0.000022X^2+
-0.046890X^1+
-16.533344
Press any key to continue . . .

Это сообщение отредактировал(а) Shanks - 7.4.2011, 12:26

ksili

Дата 7.4.2011, 12:57 (ссылка)

(нет голосов)

Загрузка ...

Эксперт

Профиль
Группа: Завсегдатай
Сообщений: 2069
Регистрация: 3.11.2005
Где: Красноярск

Репутация: 2
Всего: 17

Отклонение от чего считаете? Если от исходных точек, то значения между ними и не должны быть близкими к ним. Полиномиальную интерполяцию при большом количестве точек потому и не любят.

--------------------

Ничто так не развивает аналитическое мышление, как отладка сложной программы без возможности пошагового выполнения (с)

Shanks

Дата 7.4.2011, 13:00 (ссылка)

(нет голосов)

Загрузка ...

Новичок

Профиль
Группа: Участник
Сообщений: 4
Регистрация: 7.4.2011

Репутация: нет
Всего: нет

float x[N]={1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 3e9};
float y[N]={-17, -21, -43, -75, -96, -120, -142, -174};

Вот исходный массив. Однако в алгоритме уже при F(5500 или 5,5e3)=84.443520
дальше ещё хуже. Как заставить алгоритм работать или вообще выбрать что-то другое ?

ksili

Дата 7.4.2011, 13:13 (ссылка)

(нет голосов)

Загрузка ...

Эксперт

Профиль
Группа: Завсегдатай
Сообщений: 2069
Регистрация: 3.11.2005
Где: Красноярск

Репутация: 2
Всего: 17

Алгоритм работает. Ты скажи, чего тебе надо от него? Подозреваю, что тебе нужно не то, что делает данный алгоритм.

--------------------

Shanks

Дата 7.4.2011, 13:27 (ссылка)

(нет голосов)

Загрузка ...

Новичок

Профиль
Группа: Участник
Сообщений: 4
Регистрация: 7.4.2011

Репутация: нет
Всего: нет

Мне надо строить график примерно такого вида,а так же уметь получить значение y при любом x в этом диапазоне.
Совершенно не отрицаю, что могу пытаться использовать не тот алгоритм. Собственно для этого я и обратился за помощью к более понимающим людям.

Это сообщение отредактировал(а) Shanks - 7.4.2011, 13:56

Присоединённый файл ( Кол-во скачиваний: 27 )

gr.jpg 33,77 Kb

ksili

Дата 8.4.2011, 03:42 (ссылка)

(нет голосов)

Загрузка ...

Эксперт

Профиль
Группа: Завсегдатай
Сообщений: 2069
Регистрация: 3.11.2005
Где: Красноярск

Репутация: 2
Всего: 17

По рисунку выглядит будто там используется линейная интерполяция. Точки просто соединены отрезками прямых.
Уравнение внутри каждого отрезка от точки (x0, f(x0)) до точки (x1, f(x1)):

Цитата

y(x) = f(x0) + (x-x0)*(f(x1)-f(x0))/(x1-x0)

--------------------

Shanks

Дата 8.4.2011, 08:41 (ссылка)

(нет голосов)

Загрузка ...

Новичок

Профиль
Группа: Участник
Сообщений: 4
Регистрация: 7.4.2011

Репутация: нет
Всего: нет

Как выяснилось надо конвертить значения по оси x в логарифмические.
В общем потребовалось добавить всего-то пару строк:

Код


    for (z=0;z<N-1;z++)
    {
        x[z]=log10(x[z]);
    }

0 Пользователей читают эту тему (0 Гостей и 0 Скрытых Пользователей)
0 Пользователей:
« Предыдущая тема \| Алгоритмы \| Следующая тема »