БлогNot. Полином Ньютона на Паскале и C++

Полином Ньютона на Паскале и C++

Предупреждение: полиномы Ньютона (их два) - вещь гнилая. 1-й работает плохо в конце отрезка, 2-й - в начале. На практике вряд ли используются

Письмо:

Добрый день! Я являюсь постоянным читателем http://nickolay.info. У меня к вам огромная просьба, не могли бы вы помочь мне разобраться с задачей на Pascal связанная с численными методами, читал у вас "Реализация основных методов интерполяции функции одной переменной", но это не совсем то что мне нужно. Если бы вы помогли разобраться со следующей задачей я был бы вам очень признателен: Построить интерполяционный многочлен Ньютона по известным в узлах х0,х1,...,х5 значениями f(x). Вычислить значения этого многочлена в равномерной сетке точек [a,b] с шагом h и сравнить их со значением функции f(x) в этих точках. Использовать схему интерполяции вперед.
Функция f(x)=x+exp(-x^2).
Значения функции:
Х0=0,2
Х1=0,5
Х2=1
Х3=1,3
Х4=1,5
Х5=1.9
Отрезок:
a=0
b=2
Шаг: h=0,05

Ответ: (вдруг кому пригодится, всё ж минуты 4 "труда" :)

Добрый день, насколько я ориентируюсь на своём сайте, полином Ньютона на Паскале есть вот тут: http://nickolay.info/study/pas_methods.html оттуда и можно надёргать кода, как-то вот так, не проверял толком:

Uses Crt;
Type vector=array[1..50] of real;

Procedure Interpol (var X,Y:vector;n,nom:integer;a:real;var f,r:real);
Var
 D1,D2,D3,D4:vector;
 h,q:real;
 k,i:integer;
 Begin
      for i:=1 to n-1 do d1[i]:=y[i+1]-y[i];
      for i:=1 to n-2 do d2[i]:=d1[i+1]-d1[i];
      for i:=1 to n-3 do d3[i]:=d2[i+1]-d2[i];
      for i:=1 to n-4 do d4[i]:=d3[i+1]-d3[i];
      h:=abs(x[2]-x[1]);
      case nom of
        1:begin
            i:=1;
            while a>x[i] do i:=i+1;
            k:=i-1;
            q:=(a-x[k])/h;
            f:=y[k]+d1[k]*q+d2[k]*q*(q-1)/2+
                       d3[k]*q*(q-1)*(q-2)/6;
            R:=abs(d4[k]*q*(q-1)*(q-2)*(q-3)/24);
          end;
        2:begin
            i:=1;
            while a>x[i] do i:=i+1;
            k:=i;
            q:=(a-x[k])/h;
            f:=y[k]+d1[k-1]*q+d2[k-2]*q*(q+1)/2+
                         d3[k-3]*q*(q+1)*(q+2)/6;
            R:=abs(d4[k-4]*q*(q+1)*(q+2)*(q+3)/24);
          end;
      end;
 End;

function f(x:real):real;
begin
 f:=x+exp(-sqr(x));
end;


Var
  X,Y:vector;
  n,nom,i:integer;
  a,b,h,xx,ff,R:real;
Begin
  ClrScr;
  write('Введите размер таблицы n=');
  readln(n);
  writeln('Введите массив X=');
  for i:=1 to n do read(x[i]);
  write('массив Y=');
  for i:=1 to n do begin
   y[i]:=f(x[i]);
   write (y[i]:18:6);
  end;
  writeln;
  write('Введите интервал [a,b]=');
  readln(a,b);
  write('Введите шаг h=');
  readln(h);
  write('Введите номер интерполяционной ',
                         'формулы Ньютона: nom=');
  readln(nom);
  xx:=a;
 writeln ('X':18,'Вычислено':18,'Точное':18,'остаток':18);
 while xx<=b do begin
  Interpol(X,Y,n,nom,xx,ff,R);
  writeln(xx:18:6,ff:18:6,f(xx):18:6,R:18:6);
  xx:=xx+h;
  end;
  reset (input); readln;
End.

Скриншот тестового прогона на этих данных (Турбо Паскаль 7.1 в консольке):

Скриншот Полином Ньютона на Паскале
Скриншот Полином Ньютона на Паскале

Узлы интерполяции по оси X предполагаются равноотстоящими. Ясно, что вне интервала ]Xmin,Xmax[ вычисленные значения смысла не имеют. Также может и глюк где.

Вот C++-реализация первой интерполяционный формулы Ньютона (для левой половины интервала), запускалась в консоли Visual Studio.

#include <iostream>
using namespace std;

struct Data { double x, y; };

double u_cal (double u, int n) {
 double temp = u;
 for (int i = 1; i < n; i++) temp *= (u - i);
 return temp;
}

double fact(int n) {
 double f = 1.;
 for (int i = 2; i <= n; i++) f *= i;
 return f;
}

int main() {
 //Количество исходных точек:
 const int n = 4;
 //Исходные значения (x,y):
 Data f[n] = { 
  { 0, 2 },
  { 1, 3 },
  { 2, 12 },
  { 5, 147 } 
 }; 
 //Таблица конечных разностей (вперёд)
 double y[n][n];
 for (int i = 0; i < n; i++) y[i][0]=f[i].y;
 for (int i = 1; i < n; i++) {
  for (int j = 0; j < n - i; j++) y[j][i] = y[j + 1][i - 1] - y[j][i - 1];
 }
 //Вывод таблицы значений от x1 до x2 с шагом dx:
 double x1 = 0,
        x2 = 5,
        dx = 0.125;
 cout.width(20); cout << "X"; cout.width(20); cout << "Y" << endl;
 for (double x = x1; x <= x2; x += dx) {
  double sum = y[0][0];
  double u = (x - f[0].x) / (f[1].x - f[0].x);
  for (int i = 1; i < n; i++) {
   sum += (u_cal(u, i) * y[0][i]) / fact(i);
  }
  cout.width(20); cout.precision(14); cout << x;
  cout.width(20); cout.precision(14); cout << sum << endl;
 }

 cin.get();
 return 0;
}

Вторая интерполяционная формула Ньютона (для правой половины интервала).

#include <iostream>
using namespace std;

struct Data { double x, y; };

double u_cal (double u, int n) {
 double temp = u;
 for (int i = 1; i < n; i++) temp *= (u + i);
 return temp;
}

double fact(int n) {
 double f = 1.;
 for (int i = 2; i <= n; i++) f *= i;
 return f;
}

int main() {
 //Количество исходных точек:
 const int n = 4;
 //Исходные значения (x,y):
 Data f[n] = { 
  { 0, 2 },
  { 1, 3 },
  { 2, 12 },
  { 5, 147 } 
 }; 
 //Таблица конечных разностей (назад)
 double y[n][n];
 for (int i = 0; i < n; i++) y[i][0]=f[i].y;
 for (int i = 1; i < n; i++) {
  for (int j = n - 1; j >= i; j--) y[j][i] = y[j][i - 1] - y[j - 1][i - 1];
 }
 //Вывод таблицы значений от x1 до x2 с шагом dx:
 double x1 = 0,
        x2 = 5,
        dx = 0.125;
 cout.width(20); cout << "X"; cout.width(20); cout << "Y" << endl;
 for (double x = x1; x <= x2; x += dx) {
  double sum = y[n-1][0];
  double u = (x - f[n-1].x) / (f[1].x - f[0].x);
  for (int i = 1; i < n; i++) {
   sum += (u_cal(u, i) * y[n-1][i]) / fact(i);
  }
  cout.width(20); cout.precision(14); cout << x;
  cout.width(20); cout.precision(14); cout << sum << endl;
 }

 cin.get();
 return 0;
}

19.03.2013, 16:53 [32783 просмотра]


теги: учебное c++ алгоритм pascal

показать комментарии (1)