4. Решение функциональных и вычислительных задач химической технологии

4.1. Обработка экспериментальных данных

4.1.1 Интерполяционный многочлен Лагранжа

Интерполяционная формула Лагранжа используется для произвольно заданных узлов интерполирования.

Пусть в точках x₀, x₁, …, x_n таких, что a x₀<…< x_nb, известны значения функции y=f(x), т.е. на отрезке [a, b] задана табличная (сеточная) функция

x	x0	x1	…	xn
y	y0	y1	…	yn

Требуется построить некоторую функцию P_n(x) степени не выше n, имеющую в заданных узлах x₀, x₁, …, x_n те же значения, что и функция f(x). С этой целью представим P_n(x) в следующем виде:

где С_i – некоторая константа, которая определяется из следующего уравнения:

.

В этом случае для многочлена:

Подставляя выражение для константы С_i в исходную зависимость полинома для каждого i=0,1,2,…, n следующее значение многочлена:

,

Подставляя значение P_n(x) для всех точек от 0 до n, получим следующее выражение для многочлена Лагранжа:

Пример 4.1.1. Построить интерполяционный полином для функции y=sin x.

Возьмем сетку, состоящую из трех точек: x₀=0; x₁= ; x₂= , выпишем соответствующие этим аргументам значения функции sin x: y₀=0; y₁= ; y₂=1.

Построим по этой таблице интерполяционный полином второй степени, использую формулу Лагранжа:

.

Проверяем, что в точках сетки этот полином принимает нужные значения. Оценим погрешность интерполирования, сравнив значения sin x и интерполяционного полинома в точке х= .

.

Значительная величина погрешности определяется тем, что на отрезке длиной мы взяли грубую сетку, состоящую всего из трех точек. Чтобы улучшить точность интерполирования, нужно либо увеличить число точек n и повысить соответственно степень интерполяционного полинома P_n(x), либо уменьшить длину исходного отрезка.

4.1.2. Интерполяционный многочлен Ньютона

Предположим дополнительно, что рассматриваемые значения аргумента являются равноотстоящими, т. е. образуют арифметическую прогрессию.

В этом случае шаг таблицы h = х_i+1 - x_i (i = 0, 1, 2, ..., n) = const является величиной постоянной. Для таких таблиц построение интерполяционных формул (как впрочем, и вычисление по этим формулам) заметно упрощается.

Прежде чем перейти к рассмотрению этого вопроса, познакомимся с понятием конечных разностей.

Пусть функция задана таблицей с постоянным шагом. Разности между значениями функции в соседних узлах интерполяции называют конечными разностями первого порядка:

Δ y_i = y_i+1 - y_i (i = 0, 1, 2, ...).

Из конечных разностей первого порядка образуются конечные разности второго порядка:

Δ ²y_i = Δ yⁱ⁺¹ - Δ y_i (i = 0, 1, 2, ...)

Для записи конечных разностей используются горизонтальные и диагональные таблицы.

Горизонтальная таблица

x	y	Δy	Δ2y	Δ3y
x0	y0	Δy0	Δ2 y0	Δ3y0
x1	y1	Δy1	Δ2 y1
x2	y2	Δy2
x3	y3

Диагональная таблица

x	y	Δy	Δ2y	Δ3y
x0	y0	Δy0	Δ2 y0	Δ3y0
x1	y1
		Δy1
x2	y2		Δ2 y1
		Δy2
x3	y3

Формула Ньютона «вперед»:

где n – порядок полинома, h – шаг (расстояние между узлами).

Первый порядок:

x1=x0-h	Δy= y1 – y0
x2=x0-2h	Δy1= y2 – y1
xn=x0-n×h	Δyn-1= yn – yn-1

Второй порядок:

Δ2y= Δy1 – Δy0

Δ2y1= Δy2 – Δy1

Δ2yn-2= Δyn-1 – Δyn-2

Формула Ньютона «назад»:

Программа для интерполирования «вперед»:

Программа для интерполирования «назад»:

program newton; const n=3; var d:array[0..n,0..n] of real; x,y:array[0..n] of real; p:array[1..n] of real; f1,f2:text; h,x1,s1:real;i,j:integer; begin assign(f1,'n1.pas');assign(f2,'nnn.pas'); reset(f1);rewrite(f2);x1:=1.5;h:=1; for i:=0 to n do read(f1,x[i],y[i]); s1:=y[0]; for i:=0 to n do d[i,0]:=y[i]; for j:=1 to n do for i:=0 to n-j do d[i,j]:=d[i+1,j-1]-d[i,j-1]; for i:=1 to n do begin p[i]:=1;for j:=1 to i do p[i]:=p[i]*(x1-x[j-1])/(h*j); s1:=s1+d[0,i]*p[i]; end; write(f2,s1:7:4);close(f2);end.

program newton; const n=3; var d:array[0..n,0..n] of real; x,y:array[0..n] of real; p:array[1..n] of real; f1,f2:text; h,x1,s1:real;i,j:integer; begin assign(f1,'n1.pas');assign(f2,'nn3.pas'); reset(f1);rewrite(f2);x1:=1.5;h:=1; for i:=0 to n do read(f1,x[i],y[i]);s1:=y[n]; for i:=0 to n do d[i,0]:=y[i]; for j:=1 to n do for i:=0 to n-j do d[i,j]:=d[i+1,j-1]-d[i,j-1]; for i:=1 to n do begin p[i]:=1;for j:=1 to i do p[i]:=p[i]*(x1-x[n-j+1])/(h*j); s1:=s1+d[n-i,i]*p[i]; end; write(f2,s1:5:2);close(f2);end.