Реализации алгоритмов/Интерполяция/Многочлен Лагранжа

C#

using System;
 
namespace ConsoleApp
{
    class Program
    {
        static void Main()
        {
            const int size = 10;
            var xValues = new double[size];
            var yValues = new double[size];
 
            for (int i = 0; i < size; i++)
            {
                xValues[i] = i; 
                yValues[i] = TestF(i);
            }
 
            Console.WriteLine(InterpolateLagrangePolynomial(13.6, xValues, yValues, size));
        }
 
        static double InterpolateLagrangePolynomial (double x, double[] xValues, double[] yValues, int size)
        {
            double lagrangePol = 0;
 
            for (int i = 0; i < size; i++)
            {
                    double basicsPol = 1;
                    for (int j = 0; j < size; j++)
                    {
                        if (j != i)
                        {
                            basicsPol *= (x - xValues[j])/(xValues[i] - xValues[j]);
                        }
                    }
                    lagrangePol += basicsPol * yValues[i];
            }
 
            return lagrangePol;
        }
 
        static double TestF(double x)
        {
            return  x*x*x + 3*x*x + 3*x + 1; // for example
        }
    }
}

C

#include <stdio.h>

double InterpolateLagrangePolynomial (double x, double* x_values, double* y_values, int size)
{
	double lagrange_pol = 0;
	double basics_pol;

	for (int i = 0; i < size; i++)
	{
		basics_pol = 1;
		for (int j = 0; j < size; j++)
		{
			if (j == i) continue;
			basics_pol *= (x - x_values[j])/(x_values[i] - x_values[j]);		
		}
		lagrange_pol += basics_pol*y_values[i];
	}
	return lagrange_pol;
}

double testF(double x)
{
	return  x*x*x + 3*x*x + 3*x + 1; // for example
}

int main(void)
{
	const int size = 10;
	double x_values[size];
	double y_values[size];


	for (int i = 0; i < size; i++)
	{
		x_values[i] = i; 
		y_values[i] = testF(i);
	}

	printf ("%lf\n", InterpolateLagrangePolynomial(13.6, x_values, y_values, size));

	return 0;
}

Паскаль

var L:real;
    i,j,n:integer;
    x:array[1..10] of real;
    y:array[1..10] of real;

begin
     write('n=');readln(n);
     FOR i:=1 TO n DO
     begin
          write('x[',i,']=');readln(x[i]);
          write('y[',i,']=');readln(y[i]);
     end;
     begin
          write('x[',n+1,']=');readln(x[n+1]);
     end;
     y[n+1]:=0;
     FOR j:=1 TO n DO begin
     L:=1;
     FOR i:=1 TO n DO
     begin
         IF i<>j THEN
         begin
         L:=L*(x[n+1]-x[i])/(x[j]-x[i]);
         end;
     end;
     y[n+1]:=y[n+1]+y[J]*L;end;
     writeln('y[',n+1,']=',y[n+1]:1:0);
     FOR i:=1 TO n DO
     begin
          writeln('x[',i,']=',x[i]:10:10,' y[',i,']=',y[i]:10:10);
     end;
     begin
          writeln('x[',n+1,']=',x[n+1]:10:10,' y[',n+1,']=',y[n+1]:10:10);
     end;
     readln;
end.

Бэйсик - QBasic

DIM x(10), s(10)
INPUT "n=", n
FOR i = 1 TO n
INPUT "x(i)=", x(i)
INPUT "y(i)=", y(i)
NEXT i
INPUT "x(n+1)=", x(n + 1)

F = 0
y(n + 1) = 0
FOR j = 1 TO n
L = 1
FOR i = 1 TO n
IF i <> j THEN
        L = L * (x(n + 1) - x(i)) / (x(j) - x(i))
END IF
NEXT i
     y(n + 1) = y(n + 1) + L * y(j)
NEXT j

FOR i = 1 TO n
PRINT "x("; i; ")="; x(i); " y("; i; ")="; y(i)
NEXT i
PRINT "x("; n + 1; ")="; x(n + 1); " y("; n + 1; ")="; y(n + 1)
END

C++

double lagrange1(double* x, double* y, short n, double _x)
{
	double result = 0.0;

	for (short i = 0; i < n; i++)
	{
		double P = 1.0;

		for (short j = 0; j < n; j++)
			if (j != i)
				P *= (_x - x[j])/ (x[i] - x[j]);

		result += P * y[i];
	}	

	return result;
}

double lagrange2(double* x, double* y, short n, double _x)
{
	double result = 0.0;

	double h = x[1] - x[0];
	
	for (short i = 0; i < n; i++)
	{
		double P = 1.0;

		for (short j = 0; j < n; j++)
			if (i != j)
				P *= (_x - x[0] - h* j)/ h/ (i - j);

		result += P* y[i];
	}

	return result;
}

Octave

%Функция для которой строится многочлен.
function F = aFunc(x)
	for i=1:length(x)
		F(i) = 3.5*e^(x(i)/2) + 3.5*3^(-x(i)); %Заменить на свою, или вернуть табличное значение в данной точке.
	end;
end;

%Собственно нахождение значения интерполяционного многочлена в заданной точке. x - Точка, presetX - точки в которых известно значение функции
function R = lagrang(x, presetX)
n = length(presetX);
R=0;
for i=1:n
	li=1;
	for j=1:n
		if i ~=j
			li = li*((x-presetX(j))/(presetX(i) - presetX(j)));
		end
	end
	R=R+aFunc(presetX(i))*li;
end;

end;

Maple

lx --- список $x_{i}$ ly --- список $y_{i}(x_{i})$. n --- степень полинома. z --- здесь будет записан полином. x --- переменная в полиноме. d --- вспомогательная переменная для накопления произведения. Всё остальное смотрим в документации по Maple.

ly := [2, 6, 14]: 
lx := [0, 1, 2]:
n := 3; z := 0: d := 1: x := 'x':
for i to n do for j to n do 
    if j <> i then 
        d := (x-op(j, lx))*d/(op(i, lx)-op(j, lx)); 
    end if end do; 
    z := z+op(i, ly)*d; 
    d := 1; 
end do:

expand(z);

Ruby

def lagrange(x, x_, y_)
  x_.size.times.map do |i|
    x_.reduce(1.0) do |res, x_j|
      x_[i] == x_j ? res : res * (x - x_j) / (x_[i] - x_j)
    end * y_[i]
  end.reduce(&:+)
end

puts lagrange(5­, [-1, 0, 1], [1, 0, 1])