Метод неопределенных коэффициентов построения интерполяционного полинома
| Наиболее очевидным алгоритмом построения многочлена P(x) является нахождение его коэффициентов ai из системы уравнений | ||
 , i = 0,1,...,n
| (4) | |
| Поскольку матрица этой системы уравнений | ||
| (5) | |
| является матрицей Вандермонда и её определитель, при не совпадающих xi, не равен нулю, то система (4) имеет решение и оно единственно. Отметим, что вычисление коэффициентов полинома посредством решения системы (4) в вычислительной практике используется крайне редко. Причиной этого является плохая обусловленность матрицы (5), приводящая к заметному росту погрешности в выполнении условий интерполирования уже при сравнительно невысоких порядках полинома. К этому следует добавить, что вычислительные затраты реализации метода пропорциональны n3. Алгоритмы Ньютона и Лагранжа более устойчивы к ошибкам округления и менее трудоемки. |
Интерполяционный полином в форме Ньютона
| Интерполяционный многочлен легко определяется если его построить в виде: | ||
| Pn(x) = С0 + С1(x - x0) + C2(x - x0) (x - x1) +...+ Cn(x - x0)(x - x1)... (x - xn-1) | (5) | |
| Исходя из условия интерполяции (1) для коэффициентов Ci получим систему уравнений треугольного вида | ||
| f(x0) = С0 f(x1) = С0 + С1(x1 - x0 ) f(x2) = С0 + С1(x1 - x0 ) + C2(x2 - x0 )(x2 - x1 ) .................................................... f(xn ) = С0 + С1(xn - x0) + C2(xn - x0)(xn - x1) +...+ Cn(xn - x0)(xn - x1)... (xn- xn-1) | ||
| Из этой системы легко находятся: | ||
   и так далее.
| ||
| Величины, стоящие в правой части приведённых выше равенств, получили название разделённых разностей, соответственно, нулевого, первого и второго порядков. Для них приняты обозначения f[xi], f[xi ,xi-1], f[xi ,xi-1 ,xi-2] и т.д. С учётом этих обозначений выражение (5) можно переписать в виде: | ||
| Pn(x) = f[x0] + f[x1 ,x0](x - x0) + f[x2 ,x1 ,x0](x - x0)(x - x1) +... + f[xn ,xn-1 ,...x0](x - x0)(x - x1)...(x - xn-1) | (6) | |
| Можно показать, что | ||
| (7) | |
| Выражения (6) и (7) определяют интерполяционный полином в форме Ньютона. Вычисление полинома в Ньютоновской форме удобно при последовательном дополнении сетки (n+2) -м узлом и наращивании порядка интерполяционного полинома. При этом необходимо вычислить лишь одно дополнительное слагаемое f[xn+1 ,xn,...x0](x - x0)(x - x1)...(x - xn) в выражении (6). |
Интерполяционный полином в форме Лагранжа
| Интерполяционный полином, очевидно, можно построить в форме | ||
| (8) | |
| где gi(x) - многочлен n-ой степени, обладающий следующим свойством: | ||
 для (i,j) = 0,1,...,n, 
| (9) | |
| Свойством (8),в частности, обладает полином вида: | ||
| (10) | |
| Тогда: | 
| (11) |
где 
| (12) |
Выражения (11) и (12) совместно образуют интерполяционную формулу Лагранжа. Отметим, что коэффициенты Лагранжа (9) не зависят от значений интерполируемой функции в узлах. Это существенно снижает вычислительные затраты при интерполировании системы функций на общей сетке.
Результаты и графики
Находим r, кДж/кг методом линейной интерполяции
| t, °C | r, кДж/кг |
| ? | |
| r(t=32°C) |
Находим r, кДж/кг методом полиномиальной интерполяции
| r(x0,x1) | -3,6 | r(t=32) | |
| r(x1,x2) | -3,9 | 1137,24 | |
| r(x0,x1,x2) | -0,015 | ||
| r(x2,x3) | -4,2 | ||
| r(x1,x2,x3) | -0,015 | ||
| r(x0,x1,x2,x3) | -7,51714E-19 | ||
| r(x3,x4) | -4,5 | ||
| r(x2,x3,x4) | -0,015 | ||
| r(x1,x2,x3,x4) | -7,51714E-19 | ||
| r(x0,x1,x2,x3,x4) |
| e% |
| 0,02% |
Относительная погрешность
| Схема Горнера | |||||||
| t | r | ||||||
| 38,7 | 1107,02 | ||||||
| Пример Решения Схема Горнера | |||||||
| f(x)=1263+(38,7-0)*(-3,6+(38,7-10)*(-0,015+(38,7-20)*(-7,51714E-19+(38,7-30)*0))) |
Поэтому в качестве приближенного значения 
естественно взять интегральную сумму 
,т.е. положить:
а также
т.е 
(1)
и 
(1')
Эти приближенные равенства называются формулами прямоугольников.
В том случае, когда f(x) 
0, формулы (1) и (1’) с геометрической точки зрения означают, что площадь криволинейной трапеции aABb, ограниченной дугой кривой y=f(x), осью Ох и прямыми х=а и х=b, принимается приближенно равной площади ступенчатой фигуры, образованной из n прямоугольников с основаниями 
и высотами: y0, y1, y2, …, yn-1 – в случае формулы (1) (рис.8) и y1, y2, y3, …, yn – в случае формулы (1') (рис.9).
Исходя из приведенного выше геометрического смысла формул (1) и (1') способ приближенного вычисления определенного интеграла по этим формулам принято называть методом прямоугольников.
Всякое приближенное вычисление имеет определенную ценность лишь тогда, когда оно сопровождается оценкой допущенной при этом погрешности. Поэтому формулы прямоугольников будут практически пригодны для приближенного вычисления интегралов лишь в том случае, если будет существовать удобный способ оценки получающейся при этом погрешности (при заданном n), позволяющий к тому же находить и число частей n разбиения сегмента, гарантирующее требуемую степень точности приближенного вычисления.
Будем предполагать, что функция f(x) имеет ограниченную производную на сегменте [a, b], так что существует такое число М>0, что для всех значений х из [a, b] выполняется неравенство |f'(x)| 
M. Качественный смысл этого неравенства заключается в том, что скорость изменения значения функции ограничена. В реальных природных системах это требование практически всегда выполнено. В этих условиях абсолютная величина погрешности Rn, которую мы допускаем, вычисляя интеграл по формуле прямоугольников может быть оценена по формуле [27]:
|Rn| M(b-a)2/2n (2)
При неограниченном возрастании n выражение M(b-a)2/2n, а, следовательно, и абсолютная величина погрешности Rn будет стремиться к нулю, т.е. точность приближения будет тем больше, чем на большее число равных частей будет разделен сегмент [a, b]. Абсолютная погрешность результата будет заведомо меньше заданного числа 
>0, если взять
n > M(b-a)2/2 .
Следовательно, для вычисления интеграла с указанной степенью точности достаточно сегмент [a, b] разбить на число частей, большее числа M(b-a)2/2 .
Метод прямоугольников – это наиболее простой и вместе с тем наиболее грубый метод приближенного интегрирования. Заметно меньшую погрешность дает другой метод – метод трапеций.
Метод трапеций
Очевидно, что чем больше будет число n отрезков разбиения, тем более точный результат дадут формулы (3а) и (3б). Однако увеличение числа отрезков разбиения промежутка интегрирования не всегда возможно. Поэтому большой интерес представляют формулы, дающие более точные результаты при том же числе точек разбиения.
Простейшая из таких формул получается как среднее арифметическое правых частей формул (1) и (1'):
(4)
Легко усмотреть геометрический смысл этой формулы. Если на каждом отрезке разбиения дугу графика подинтегральной функции y=f(x) заменить стягивающей ее хордой (линейная интерполяция), то мы получим трапецию, площадь которой равна 
и следовательно, формула (4) представляет собой площадь фигуры, состоящей из таких трапеций (рис.10). Из геометрических соображений понятно, что площадь такой фигуры будет, вообще говоря, более точно выражать площадь криволинейной трапеции, нежели площадь ступенчатой фигуры, рассматриваемая в методе прямоугольников.
Приведя в формуле (4) подобные члены, окончательно получим
(5)
Формулу (5) называют формулой трапеций.
Формулой трапеций часто пользуются для практических вычислений. Что касается оценки погрешности Rn, возникающей при замене левой части (5) правой, то доказывается, что абсолютная величина ее удовлетворяет неравенству:
(6)
где М2 – максимум модуля второй производной подинтегральной функции на отрезке [a,b], т.е.
.
Следовательно, Rn убывает при 
по крайней мере так же быстро, как 
.
Абсолютная погрешность Rn будет меньше наперед заданного числа > 0, если взять 
.
Метод парабол (метод Симпсона)
Значительное повышение точности приближенных формул может быть достигнуто за счет повышения порядка интерполяции. Одним из таких методов приближенного интегрирования является метод парабол. Идея метода исходит из того, что на частичном промежутке дуга некоторой параболы в общем случае теснее прилегает к кривой y=f(x), чем хорда, соединяющая концы дуги этой кривой, и поэтому значения площадей соответствующих элементарных трапеций, ограниченных “сверху” дугами парабол, являются более близкими к значениям площадей, соответствующих частичных криволинейных трапеций, ограниченных сверху дугой кривой y=f(x), чем значения площадей соответствующих прямолинейных трапеций. Сущность метода заключается в следующем. Отрезок [a,b] делится на 2n равных частей. Пусть точки деления будут
х0=а, x1, x2, …x2n-2, x2n-1, x2n=b,
а y0, y1, …y2n – соответствующие значения подинтегральной функции на отрезке [a,b]. Произведем квадратичную интерполяцию данной подинтегральной функции на каждом из отрезков разбиения (заменим дугу графика подинтегральной функции дугой параболы с вертикальной осью) (рис.11).
Приведем без вывода формулу парабол в окончательном виде:
(7)
(Подробный вывод формулы (7) см. в [13]).
Если подинтегральная функция f(x) имеет на отрезке [a,b] непрерывную четвертую производную, то для поправочного члена формулы (7) имеет место оценка
(8)
где М4 - максимум модуля четвертой производной подинтегральной функции на отрезке [a,b].
Сравнивая между собой оценки (6) и (8), замечаем, что с увеличением n поправочный член формулы трапеций уменьшается пропорционально величине , а для формулы парабол – пропорционально величине 
, т.е. метод парабол сходится значительно быстрее метода трапеций, тогда как с точки зрения техники вычислений оба метода одинаковы.
Результаты и график
Находим точное значение интеграла, используя первообразную.
Находим численное значение интеграла для n = 10, используя метод трапеций и метод Симпсона. Программная реализация в Pascal.
program rrr;
uses crt;
var a,b,h,x,y,n,I,I01,I02,I0,y1,I1:real;
k,l:integer;
function f1 (x:real):real;
Begin
f1:=sqrt(x*x-36);
end;
Begin
clrscr;
write('vvedite nizhnyiu granicy a=');
readln(a);
write('vvedite verhnyiu granicy b=');
readln(b);
write('vvedite wag h=');
readln(h);
writeln('kol-vo wagov posle zapiatoi l=');
readln(l);
writeln(a,b,h,l);
I01:=0.5*b*sqrt(b*b-36)-18*Ln(sqrt(b*b-36)+b);
I02:=0.5*a*sqrt(a*a-36)-18*Ln(sqrt(a*a-36)+a);
I0:=I01-I02;
writeln('I0=',I0:1:l);
x:=a+h;
n:=b-h;
k:=0;
y:=0;
while (x<=n) do
begin k:=k+1;
if (k mod 2=0) then y:=y+2*f1(x) else y:=y+4*f1(x);
x:=x+h;
end;
y:=y+f1(a)+f1(b);
I:=(h/3)*(y);
writeln('I=',I:1:l);
y1:=0;
x:=a+h;
n:=b-h;
Begin
while (x<=n) do
Begin
y1:=y1+2*f1(x);
x:=x+h;
end;
y1:=y1+f1(a)+f1(b);
I1:=(h/2)*y1;
writeln('I1=',I1:1:l);
end;
readln;
end.
| x | f(x) | Вычисления в MathCad | ||||
| 6,5 | 2,5 | S | ||||
| 6,6 | 2,749545 | 8,853 | ||||
| 6,7 | 2,98161 | |||||
| 6,8 | 3,2 | Вычисления в Pascal | ||||
| 6,9 | 3,407345 | Метод Симпсона | Погрешность, % | |||
| 3,605551 | S | |||||
| 7,1 | 3,796051 | 8,853 | ||||
| 7,2 | 3,97995 | |||||
| 7,3 | 4,158125 | Метод Трапеций | Погрешность, % | |||
| 7,4 | 4,331282 | S | ||||
| 7,5 | 4,5 | 8,853 | ||||
| 7,6 | 4,664762 | |||||
| 7,7 | 4,825971 | |||||
| 7,8 | 4,983974 | |||||
| 7,9 | 5,139066 | |||||
| 5,291503 | ||||||
| 8,1 | 5,441507 | |||||
| 8,2 | 5,589275 | |||||
| 8,3 | 5,73498 | |||||
| 8,4 | 5,878775 | |||||
| 8,5 | 6,020797 |
S=((2.5+6,02)/2)*2=8,852