Интерполяционный полином Лагранжа на Паскале и C++

Для вычисления полинома Лагранжа в точке x0 по заданным значениям x[i],y[i] достаточно следующей функции:

function lagrange (x0:real; var x,y:array of real):real;
var k,j:integer; l,s:real;
begin
 l:=0;
 for k:=low(y) to high(y) do begin
  s:=1;
  for j:=low(x) to high(x) do begin
   if j<>k then s:=s*(x0-x[j])/(x[k]-x[j]);
  end;
  l:=l+y[k]*s;
 end;
 lagrange:=l;
end;

Для проверки заполним данные x[i],y[i] значениями синуса на интервале от 0 до pi и затем рассчитаем с помощью нашей функции значения на границах и серединах интервалов.

var x,y:array [1..11] of real;
 x0,y0,y1:real;
 i:integer;
begin
 x[1]:=0;
 for i:=2 to 11 do begin
  x[i]:=x[i-1]+pi/10;
  y[i]:=sin(x[i]);
 end;
 x0:=0;
 writeln('X':8,'SIN(X)':8,'L(X)':8);
 for i:=1 to 21 do begin
  y0:=sin(x0);
  y1:=lagrange(x0,x,y);
  writeln (x0:8:2,y0:8:2,y1:8:2);
  x0:=x0+pi/20;
 end;
 write ('ENTER to exit');
 readln;
end.

Так как всё совпало, полином реализован верно.

Тест:

X SIN(X) L(X)
0.00 0.00 0.00
0.16 0.16 0.16
0.31 0.31 0.31
0.47 0.45 0.45
0.63 0.59 0.59
0.79 0.71 0.71
0.94 0.81 0.81
1.10 0.89 0.89
1.26 0.95 0.95
1.41 0.99 0.99
1.57 1.00 1.00
1.73 0.99 0.99
1.88 0.95 0.95
2.04 0.89 0.89
2.20 0.81 0.81
2.36 0.71 0.71
2.51 0.59 0.59
2.67 0.45 0.45
2.83 0.31 0.31
2.98 0.16 0.16
3.14 0.00 0.00
ENTER to exit

А вот метод Лагранжа на С++, строим полином порядка (n-1) по n точкам (x,y) без обращения матрицы. Предполагается, что данные корректны. Проверялось в консоли Visual Studio 2015.

#include <iostream>
using namespace std;

struct Data { double x, y; };

double interpolate(Data f[], int n, double xi) {
 double result = 0;
 for (int i = 0; i < n; i++) {
  double term = f[i].y;
  for (int j = 0; j < n; j++) {
   if (j != i) term *= (xi - f[j].x) / (f[i].x - f[j].x);
  }
  result += term;
 }
 return result;
}

int main() {
 //Количество исходных точек:
 const int n = 4;
 //Исходные значения (x,y):
 Data f[n] = { 
  { 0, 2 },
  { 1, 3 },
  { 2, 12 },
  { 5, 147 } 
 }; 
 //Вывод таблицы значений от 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 y = interpolate(f, n, x);
  cout.width(20); cout.precision(14); cout << x;
  cout.width(20); cout.precision(14); cout << y << endl;
 }

 cin.get();
 return 0;
}

24.11.2013, 14:18 [26569 просмотров]