Главная страница Случайная страница
КАТЕГОРИИ:
АвтомобилиАстрономияБиологияГеографияДом и садДругие языкиДругоеИнформатикаИсторияКультураЛитератураЛогикаМатематикаМедицинаМеталлургияМеханикаОбразованиеОхрана трудаПедагогикаПолитикаПравоПсихологияРелигияРиторикаСоциологияСпортСтроительствоТехнологияТуризмФизикаФилософияФинансыХимияЧерчениеЭкологияЭкономикаЭлектроника
Интерполяционный полином Ньютона
|
|
Ньютон предложил интерполяционный полином степени n в виде:
P n (x)= A0 +A1 (x-x 0 ) +A2 (x-x 0 )(x-x 1 ) +...+An (x-x 0 )(x-x 1 )...(x-x n-1 ). (17)
Коэффициенты полинома Ньютона Ai определяются из условий:
P n (x i ) = f i, i=0, 1, 2...n.
Полагаем x = x 0, тогда в формуле (17) все слагаемые, кроме A0, обращаются в нуль, следовательно, A0 = f 0.
Затем полагаем x = x 1, тогда f 0 + A1 (x 1 –x 0 ) = f 1, откуда находим коэффициент
A1 = (f 0- f 1 ) / (x 0 -x 1 ) = f 01,
который называется разделенной разностью первого порядка.
При x = x2 полином принимает значение
P n (x 2 )=f 0 + f 01 (x 2 -x 0 ) + A2 (x 2 -x 0 )(x 2 -x 1 ),
поэтому
A2 = (f 01- f 02 ) / (x 1 -x 2 ) = f 012, где f 02 = (f 0- f 2 ) / (x 0 -x 2 ).
Величина f 012 называется разделенной разностью второго порядка. Аналогичным образом при x = x3 находим коэффициент полинома A3 и т.д. Для коэффициента Ak имеем следующее выражение:
Ak = (f 01…k-1- f 01…k ) / (x k-1 -x k ).
Результаты расчета разделенных разностей сведем в таблицу:
| x | f(x) | ||||
| x 0 | f 0 | ||||
| x 1 | f 1 | f 01 = (f 0- f 1 ) / (x 0 -x 1 ) | |||
| x 2 | f 2 | f 02 = (f 0- f 2 ) / (x 0 -x 2 ) | f012 = (f 01- f 02 ) / (x 1 -x 2 ) | ||
| x 3 | f 3 | f 03 = (f 0- f 3 ) / (x 0 -x 3 ) | f 013 = (f 01- f 03 ) / (x 1 -x 3 ) | f 0123 = (f 012- f 013 ) / (x 2 -x 3 ) | |
| x 4 | f 4 | f 04 = (f 0- f 4 ) / (x 0 -x 4 ) | f 014 = (f 01- f 04 ) / (x 1 -x 4 ) | f 0124 = (f 012- f 014 ) / (x 2 -x 4 ) | f 01234 = (f 0123- f 0124 ) / (x 3 -x 4 ) |
Для построения интерполяционного полинома Ньютона используются только диагональные элементы этой таблицы, остальные элементы являются промежуточными данными.
После определения коэффициентов полинома Ньютона вычисление его значений при конкретных аргументах x наиболее экономично проводить по схеме Горнера, получаемой путем последовательного вынесения за скобки множителей (x-x i ) в формуле (17):
P n (x) = f 0 +(x-x 0 )(f 01 +(x-x 1 )(f 012 +(x-x 2 )(f 0123 +...)...). (18)
В отличие от алгоритма вычисления полинома Лагранжа при интерполяции полиномом Ньютона удается разделить задачи определения коэффициентов и вычисления значений полинома при различных значениях аргумента, что ускоряет вычисления.
Ниже приведена программа интерполяции полиномом Ньютона. Вычисление коэффициентов полинома Ньютона осуществляется в строках 110-150 программы с помощью двух вложенных циклов по строкам (переменная I) и столбцам (переменная J) вышеприведенной таблицы.
Разделенные разности размещаются в массиве Y(X), где первоначально хранились значения функции в узлах. При вычислении разделенных разностей очередного порядка этот массив частично обновляется.
Для контроля и отладки программы использован тот же пример, что и при отладке программы интерполяции полиномом Лагранжа.
Программа интерполяции полиномом Ньютона:
10 PRINT " Интерполяция полиномом Ньютона для N+1 узлов."
20 INPUT " Введите N="; N: DIM X(N), Y(N)
30 PRINT " Введите значения X(I) и Y(I), I=0, 1, 2...N."
40 FOR I=0 TO N: PRINT " X" I; ", Y" I;: INPUT X(I), Y(I): NEXT I
50 GOSUB 110
60 PRINT: PRINT " Коэффициенты полинома Ньютона: "
70 FOR I=0 TO N: PRINT " A" I; " =" Y(I): NEXT I
80 PRINT: INPUT " Введите X="; X1
90 GOSUB 160
100 PRINT " Y(X)=" P: GOTO 80
110 REM Вычисление коэффициентов полинома Ньютона
120 FOR J=1 TO N: A=Y(J-1): B=X(J-1)
130 FOR I=J TO N: Y(I)=(A-Y(I))/(B-X(I)): NEXT I
140 NEXT J
150 RETURN
160 REM Вычисление полинома
170 P=Y(N)
180 FOR I=N-1 TO 0 STEP -1: R=X1-X(I): P=Y(I)+R*P: NEXT I
190 RETURN
Расчет на ЭВМ:
Интерполяция полиномом Ньютона для N+1 узлов.
Введите N=? 5
Введите значения X(I) и Y(I), I=0, 1, 2...N.
X 0, Y 0? 1.00, 0.84270
X 1, Y 1? 1.05, 0.86244
X 2, Y 2? 1.10, 0.88021
X 3, Y 3? 1.15, 0.89612
X 4, Y 4? 1.20, 0.91031
X 5, Y 5? 1.25, 0.92290
Коэффициенты полинома Ньютона:
A 0 = 0.84270
A 1 = 0.39480
A 2 = -0.39400
A 3 = 0.14682
A 4 = 0.19800
A 5 = -1.31426
Введите X=? 1.175
Y(X) = 0.90342
Введите X=?
Обратная интерполяция состоит в определении значения аргумента x, которому соответствует заданное значение функции y.
Задачу обратной интерполяции легко обратить, считая значения функции, наоборот, значениями аргумента.
4. Численное интегрирование
Ставится задача вычислить интеграл вида
где a и b - нижний и верхний пределы интегрирования; f(x) - непрерывная функция на отрезке [ a, b ].
К численному интегрированию обращаются, когда нельзя через элементарные функции аналитически записать первообразную интеграла или когда подобная запись имеет сложный вид.
Общий подход к решению задачи является следующим. Определенный интеграл J представляет собой площадь, ограниченную кривой f(x) и прямыми x=a и x=b. Эту площадь вычисляют, разбивая интервал от a до b на множество мелких интервалов, находя приблизительно площадь каждой полоски, получающейся при таком разбиении, и суммируя площади этих полосок.
Сущность большинства методов вычисления определенных интегралов состоит в замене подынтегральной функции f(x) аппроксимирующей функцией j(x), для которой можно легко записать первообразную в элементарных функциях, т.е.
где S - приближенное значение интеграла; R - погрешность вычисления интеграла.
Рассмотрим метод Симпсона, который получил широкое применение на практике при численном интегрировании с использованием ЭВМ (рис. 4).
Подынтегральную функцию f(x) заменим интерполяционным полиномом второй степени P 2 (x) - параболой, проходящей через узлы x 0, x 1, x 2, тогда
где R - погрешность вычисления интеграла.
Для записи полинома P 2 (x) воспользуемся интерполяционной формулой Ньютона (17) для трех узлов:
P 2 (x) = f 0 + f 01 (x - x 0 ) + f 012 (x - x 0 )(x - x 1 ), где
f 01 = (f 0 - f 1 )/(x 0 - x 1 )=(f 1 -f 0 )/h,
f 012 =(f 01 - f 02 )/(x 1 – x 2 )=(f 0 - 2 f 1 +f 2 )/ 2 h,
где f 01, f 012 - разделенные разности; h - расстояние между узлами.
Введем новую переменную z = x – x 0, тогда x = z + x 0 и полином принимает вид:
P 2 (z)=f 0 + (f 01 - f 012 h)z + f 012 z 2.
Теперь вычислим интеграл от полинома:
(19)
Из формулы (19) получим формулу для приближенного вычисления интеграла по всему интервалу [ a, b ]. Для этого разобьем интервал [ a, b ] на m интервалов длиной 2 h и для каждого из них применим формулу (19).
В результате получим:
(20)
Эту формулу называют формулой Симпсона или формулой парабол. Погрешность вычисления интеграла по формуле Симпсона определяется выражением:
где h - промежуточная точка из интервала (0, h).
Часто для оценки погрешности формулы Симпсона пользуются следующим приемом: последовательно вычисляют интеграл при делении интервала [ a, b ] на m и 2 m частей. Если получающиеся при этом значения обозначить соответственно через Sm и S 2 m, то совпадение первых знаков Sm и S 2 m позволяет судить о точности полученных значений.
Ниже приведена программа численного интегрирования методом Симпсона. В качестве примера применения метода она использована для вычисления функции ошибок, которая выражается интегралом:
После введения пределов и числа интервалов интегрирования в строке 50 вычисляется расстояние между узлами h. В строках 60-100 реализуется вычисление интеграла по формуле Симпсона. После умножения в строке 105 интеграла на коэффициент в строке 110 производится вывод его значения на экран дисплея. Значение подынтегральной функции вычисляется с помощью подпрограммы в строке 120.
Далее приведены результаты расчета на ЭВМ значения функции ошибок при аргументе x =0.678. Табличное значение функции равно erf( 0.678 ) =0.66236. Расчет показывает, что при числе интервалов интегрирования M=4 погрешность вычисления интеграла не превышает 10-5.
Программа численного интегрирования методом Симпсона:
10 PRINT " Численное интегрирование методом Симпсона."
20 INPUT " Введите нижний предел интегрирования A="; A
30 INPUT " Введите верхний предел интегрирования B="; B
40 INPUT " Введите число интервалов интегрирования M="; M
50 H=(B-A)/M/2
60 X=A: GOSUB 120: J=F: N=0
70 X=X+H: GOSUB 120: J=J+4*F
80 N=N+2: IF N=2*M THEN 100
90 X=X+H: GOSUB 120: J=J+2*F: GOTO 70
100 X=B: GOSUB 120: J=(J+F)*H/3
105 PI=3.14159265: J=(2/SQR(PI))*J
110 PRINT " Значение интеграла J="; J: GOTO 40
120 F=EXP(-(X, 2)): RETURN
130 END
Расчет на ЭВМ:
Численное интегрирование методом Симпсона.
Введите нижний предел интегрирования A=? 0
Введите верхний предел интегрирования B=? 0.678
Введите число интервалов интегрирования M=? 2
J = 0.6623778
Введите число интервалов интегрирования M=? 4
J = 0.6623601
Вычисление кратных интегралов можно осуществить повторным применением алгоритмов и программ численного интегрирования. Однако, с повышением кратности интегралов резко возрастает объем вычислительной работы. В этом случае целесообразно использовать методы Монте-Карло (методы статистических испытаний), которые свободны от этого недостатка, хотя и обеспечивают сравнительно невысокую точность. Один из вариантов метода Монте-Карло представлен на рис. 5.
Производится n испытаний, при которых берется случайное число x i, рав-номерно распределенное в интервале интегрирования [ a, b ], и вычисляется площадь прямоугольника f(x i )× (b-a), которая принимается за приближенное значение интеграла функции f(x). Вследствие случайности значения x iпогрешность интеграла также будет носить случайный характер. Приближенное значение интеграла вычисляется как среднее после n испытаний:
(21)
Погрешность интеграла будет уменьшаться с ростом числа испытаний n по закону e ~ n-1/2.
Формула (21) обобщается на кратные интегралы:
, (22)
где v s - k-мерный объем области интегрирования.
Число узлов i, в которых необходимо вычислить подынтегральную функцию, пропорционально e -2 независимо от кратности интеграла.
Ниже приведена программа численного интегрирования методом Монте-Карло. В качестве датчика псевдослучайных чисел используется стандартная функция Бейсика RND(X), генерирующая равномерно распределенные числа на интервале [0, 1]. Случайное значение аргумента задается в программе выражением A+RND(1)*H, где H=(B-A) - ширина интервала интегрирования. Так как значение функции RND(X) не зависит от значения аргумента, то он взят постоянным, равным единице.
В качестве примера программа использована для вычисления того же интеграла, что и в методе Симпсона. Погрешность порядка 10-2 обеспечивается при n > 1000.
Программа численного интегрирования методом Монте-Карло:
10 PRINT " Численное интегрирование методом Монте-Карло."
20 INPUT " Введите нижний предел интегрирования A="; A
30 INPUT " Введите верхний предел интегрирования B="; B
40 INPUT " Введите число испытаний N="; N
50 H=B-A: S=0
60 FOR I=1 TO N: S=S+EXP(-(A+RND(1)*H), 2): NEXT I
70 PI=3.14159265: S=(2/SQR(PI))*H*S/N
80 PRINT " Значение интеграла J=" S: GOTO 40
90 END
Расчет на ЭВМ:
Численное интегрирование методом Монте-Карло.
Введите нижний предел интегрирования A=? 0
Введите верхний предел интегрирования B=? 0.678
Введите число испытаний N=? 10
J = 0.6397168
Введите число испытаний N=? 100
J = 0.6517046
Введите число испытаний N=? 1000
J = 0.6634888
Введите число испытаний N=?
Рис.4. Графическая иллюстрация метода Симпсона
Рис.5. Графическая иллюстрация метода Монте-Карло
5. Интерполяция сплайнами
Интерполяция полиномами не всегда дает удовлетворительные результаты при аппроксимации зависимостей. Например, это имеет место при представлении резонансных кривых колебательных систем.
В этом случае может быть использована интерполяция сплайнами (spline – линейка, рейка). Наиболее распространена интерполяция кубическими сплайнами. При сплайновой интерполяции на каждом интервале [ x i-1, x i] строится отдельный полином третьей степени со своими коэффициентами:
j i (x) = a i + b i (x - x i-1 ) + c i (x - x i-1 ) 2+ d i (x - x i-1 ) 3.
Коэффициенты сплайнов находятся из условий «сшивания» соседних сплайнов в узловых точках:
1) равенство значений сплайнов и аппроксимируемой функции узловых точках: j i (x i-1 )= f i-1, j i (x i )= f i;
2) непрерывность первой и второй производных от сплайнов в узлах:
j’ i (x i ) = j’ i+1 (x i ), j’’ i (x i ) = j’’ i+1 (x i );
3) в соответствии с законом упругости четвертая производная от сплайнов в узлах равна нулю: j IVi (x i ).
В результате получают систему линейных алгебраических уравнений, ре-шая которую находят коэффициенты полиномов.