Эксперименты лаба4,5(2курс)
.pdf3.Аппроксимация функций
3.1.Постановка задачи
Пусть некоторая величина y является функцией аргумента x, но явная связь между y и x неизвестна (либо известная зависимость y = f(x) слишком громоздка для численных расчетов). Допустим, что в результате экспериментов получена таблица значений fxi; yig, требуется же получить значения y в других точках, отличных от узлов xi. Эта проблема решается в задаче о приближении (аппроксимации) функции: функцию f(x), явный вид которой неизвестен, требуется приближенно заменить некоторой функцией '(x) (наз. аппроксимирующей), так чтобы отклонение от f(x) в заданной области было наименьшим. Построенная таким образом аппроксимация называется точечной (примеры: интерполирование, среднеквадратичное приближение и т.д.)
Одним из основных типов точечной аппроксимации является интерполирование: для заданной функции f(x) строится интерполирующая функция '(x), принимающая в заданных точках xi те же значения yi, что и функция f(x):
'(xi) = yi; i = 0; 1; ::: ; n (1)
причем xi 6= xk при i 6= k, xi – узлы интерполяции.
Интерполирующая функция может строиться сразу для всего рассматриваемого интервала x – глобальная интерполяция, или отдельно для разных частей этого интервала – кусочная (локальная) интерполяция. Если полученная функция '(x) применяется для нахождения значения функции f(x) за пределами отрезка, содержащего узлы, то говорят об экстраполяции.
Рассмотрим использование в качестве функции '(x) интерполяцион-
ного многочлена |
|
'(x) = Pm(x) = a0 + a1x + a2x2 + ::: + amxm: |
(2) |
При глобальной интерполяции мы будем использовать все n + 1 уравнений системы (1), что позволяет найти n+1 коэффициент, откуда следует, что максимальная степень интерполяционного многочлена – m = n:
Pn(x) = a0 + a1x + a2x2 + ::: + anxn: |
(3) |
Подставляя (3) в (1) получаем:
a0 + a1x0 + a2x20::: + anxn0 = y0 a0 + a1x1 + a2x21::: + anxn1 = y1
::: |
(4) |
a0 + a1xn + a2x2n::: + anxnn = yn
(4) – система линейных алгебраических уравнений относительно неизвестных коэффициентов ai. Определитель такой системы отличен от нуля, если среди узлов xi нет совпадающих. Следовательно, в этом случае система (4) имеет единственное решение. Решив систему (4), построим интерполяционный многочлен. Такой метод построения носит название метода неопределенных коэффициентов.
Недостатки метода:
–при большом количестве узлов получается высокая степень многочлена,
–привязка к узлам интерполяции, которые, если они получены в результате измерений, могут содержать случайные погрешности.
Другой способ – подбор наиболее простой аппроксимирующей функции, график которой проходит максимально близко от узлов.
Мера отклонения функции '(x) от заданной функции f(x):
n |
|
Xi |
|
S = j'(xi) ¡ yij2: |
(5) |
=0 |
|
Метод наименьших квадратов состоит в подборе аппроксимирующей функции так, чтобы S было наименьшим.
3.2.Интерполирование
Кусочно-линейная интерполяция
КЛИ состоит в том, что заданные точки |
|
(xi; yi) соединяются прямолинейными отрезка- |
|
ми, и функция f(x) приближается ломанной с |
|
вершинами в узлах. Всего имеется n интервалов |
|
(xi¡1; xi), для каждого из них интерполяцион- |
|
ным многочленом является уравнение прямой, |
|
проходящей через две точки. |
|
Например, для i-го интервала уравнение |
|
прямой, проходящей через точки (xi¡1; yi¡1) и |
|
(xi; yi) имеет вид: |
Рис. 1. КЛИ |
|
2
|
|
y ¡ yi¡1 |
= |
x ¡ xi¡1 |
: |
|
|
|
(6) |
|||
Отсюда находим |
|
yi ¡ yi¡1 |
|
xi ¡ xi¡1 |
|
|
|
|
||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
y = aix + bi; |
|
xi¡1 6 x 6 xi; |
(7) |
||||||||
a = |
yi ¡ yi¡1 |
; |
b |
= y |
i¡1 ¡ |
a |
x |
i¡1 |
: |
|||
i |
xi ¡ xi¡1 |
|
i |
|
i |
|
|
Многочлен Лагранжа
Будем строить интерполяционный многочлен, единый для всего отрезка [x0; xn], в виде линейной комбинации многочленов степени n
L(x) = y0l0(x) + y1l1(x) + ::: + ynln(x) |
(8) |
так, чтобы многочлены li(x) обращались в нуль во всех узлах интерполяции, кроме i -го, где он должен равняться единице. Этим условиям при i = 0 отвечает многочлен вида:
l0(x) = |
(x ¡ x1)(x ¡ x2)(x ¡ x3):::(x ¡ xn) |
: |
|
(x0 ¡ x1)(x0 ¡ x2)(x0 ¡ x3):::(x0 ¡ xn) |
|||
|
|
||
Аналогично, |
|
|
|
l1(x) = |
(x ¡ x0)(x ¡ x2)(x ¡ x3):::(x ¡ xn) |
; |
|
(x1 ¡ x0)(x1 ¡ x2)(x1 ¡ x3):::(x1 ¡ xn) |
|||
|
|
||
l2(x) = |
(x ¡ x0)(x ¡ x1)(x ¡ x3):::(x ¡ xn) |
; |
|
(x2 ¡ x0)(x2 ¡ x1)(x2 ¡ x3):::(x2 ¡ xn) |
|||
::: |
|
||
|
|
(9)
(10)
li(x) = |
(x ¡ x0):::(x ¡ xi¡1)(x ¡ xi+1):::(x ¡ xn) |
|
; |
::: |
(xi ¡ x0):::(xi ¡ xi¡1)(xi ¡ xi+1):::(xi ¡ xn) |
|
|
|
|
|
|
ln(x) = |
(x ¡ x0)(x ¡ x1)(x ¡ x2):::(x ¡ xn¡1) |
: |
|
|
(xn ¡ x0)(xn ¡ x1)(xn ¡ x2):::(xn ¡ xn¡1) |
|
|
Подставляя (9),(10) в (8) получаем интерполяционный многочлен Лагран-
жа: |
n |
|
(x ¡ x0):::(x ¡ xi¡1)(x ¡ xi+1):::(x ¡ xn) |
|
|
||||||||
|
|
|
|
||||||||||
|
L(x) = |
yi |
|
: (11) |
|||||||||
|
Xi |
|
(xi |
¡ |
x0):::(xi |
¡ |
xi 1)(xi |
¡ |
xi+1):::(xi |
¡ |
xn) |
|
|
|
|
|
|
¡ |
|
|
|
|
|||||
|
=0 |
|
|
|
|
|
|
|
|
|
|
|
|
3
Единственность найденного решения следует из единственности решения системы (4). Если положить n = 1 в (11), то получим рассмотренный ранее случай линейной интерполяции, при n = 2 – случай квадратичной интерполяции
L(x) = |
(x ¡ x1)(x ¡ x2) |
y0 + |
(x ¡ x0)(x ¡ x2) |
y1+ |
|
|
(x0 ¡ x1)(x0 ¡ x2) |
(x1 ¡ x0)(x1 ¡ x2) |
|
||||
|
|
|
|
|||
|
|
|
+ |
(x ¡ x0)(x ¡ x1) |
y2: (12) |
|
|
|
|
(x2 ¡ x0)(x2 ¡ x1) |
|||
|
|
|
|
|
Многочлен Ньютона
При построении интерполяционного многочлена Лагранжа не накладывалось никакого требования на распределение узлов интерполяции. Рассмотрим случай равноотстоящих по оси x узлов интерполяции. Введем h = xi ¡ xi¡1 – шаг интерполяции, h = const.
Конечные разности
Разности первого порядка (первые разности):
¢y0 = y1 ¡ y0 = f(x0 + h) ¡ f(x0);
¢y1 = y2 ¡ y1 = f(x0 + 2h) ¡ f(x0 + h);
:::
¢yn¡1 = yn ¡ yn¡1 = f(x0 + nh) ¡ f(x0 + (n ¡ 1)h):
Разности второго порядка (вторые разности):
¢2y0 = ¢y1 ¡ ¢y0; ¢2y1 = ¢y2 ¡ ¢y1;
Разности порядка k:
¢kyi = ¢k¡1yi+1 ¡ ¢k¡1yi; i = 0; 1; :::; n ¡ 1:
Конечные разности выражаются через значения функции
¢2y0 = ¢y1 ¡ ¢y0 = (y2 ¡ y1) ¡ (y1 ¡ y0) = y2 ¡ 2y1 + y0;
¢3y0 = ¢2y1 ¡ ¢2y0 = ::: = y3 ¡ 3y2 + 3y1 ¡ y0:
В общем случае
k |
|
Xj |
|
¢kyi = (¡1)jCkj yi+k¡j; |
(13) |
=0 |
|
4
Cj |
= |
|
k! |
: |
|
|
|||
|
|
|||
k |
|
j!(k ¡ j)! |
|
|
|
|
|
Интерполяционный многочлен Ньютона будем искать в следующем виде:
N(x) = a0 + a1(x ¡ x0) + a2(x ¡ x0)(x ¡ x1) + ::: +
+ an(x ¡ x0)(x ¡ x1):::(x ¡ xn¡1): (14)
График многочлена должен проходить через заданные узлы, то есть N(xi) = yi i = 0; 1; :::; n. Для нахождения коэффициентов многочлена получаем систему
N(x0) = a0 = y0;
N(x1) = a0 + a1(x1 ¡ x0) = a0 + a1h = y1;
N(x2) = a0 + a1(x2 ¡ x0) + a2(x2 ¡ x0)(x2 ¡ x1) = = a0 + 2a1h + 2a2h2 = y2;
::: |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Отсюда находим коэффициенты ai: |
|
|
|
|
|
|
|
|
|
|
|
|
|
||||||||||||||||
|
a0 = y0; a1 = |
y1 ¡ a0 |
= |
y1 ¡ y0 |
= |
|
¢y0 |
; |
|
||||||||||||||||||||
|
|
|
|
h |
h |
|
h |
|
|||||||||||||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||||||||
a2 = |
y2 ¡ a0 ¡ 2a1h |
= |
|
y2 ¡ a0 ¡ 2¢y0 |
|
= |
|
¢2y0 |
; |
||||||||||||||||||||
|
|
|
|
|
|
|
|
|
|
|
2h2 |
||||||||||||||||||
|
|
|
|
2h2 |
|
|
|
|
|
|
|
|
::: |
|
2h2 |
|
|
|
|
|
|
|
|||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||
|
|
|
|
|
|
¢ky0 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||||
|
|
|
|
ak = |
|
|
; |
|
|
|
k = 0; 1; 2; :::; n: |
|
|
|
|
|
|
||||||||||||
|
|
|
|
k!hk |
|
|
|
|
|
|
|||||||||||||||||||
Подставляя в (14), получаем |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||||||||||
|
¢y0 |
|
|
|
|
¢2y0 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||||
N(x) = y0 + |
|
(x ¡ x0) + |
|
|
|
|
(x ¡ x0)(x ¡ x1) + :::+ |
|
|||||||||||||||||||||
h |
2!h2 |
|
|||||||||||||||||||||||||||
|
|
|
|
|
|
|
+ |
|
¢ny0 |
(x ¡ x0)(x ¡ x1):::(x ¡ xn¡1): (15) |
|||||||||||||||||||
|
|
|
|
|
|
|
|
n!hn |
|||||||||||||||||||||
Если ввести переменную t = |
x ¡ x0 |
, то |
|
|
|
|
|
|
|
|
|
||||||||||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
h |
|
|
|
|
|
|
|
|
|
|
|
|
|
||
|
x = x |
0 |
+ th; |
x ¡ x1 |
= |
x ¡ x0 ¡ h |
|
= t |
¡ |
1; |
|
||||||||||||||||||
|
|
|
|
|
|
|
h |
|
|
|
|
|
h |
|
|
|
|
|
|
|
|||||||||
|
x ¡ x2 |
= t |
¡ |
2; :::; |
|
|
x ¡ xn¡1 |
= t |
¡ |
n + 1: |
|
||||||||||||||||||
|
h |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
h |
|
|
|
|
|
|
5
В результате получаем
N(x) = N(x0 + th) = y0 + t¢y0 + t(t ¡ 1)¢2y0 + ::: + 2!
+ t(t ¡ 1):::(t ¡ n + 1)¢ny0: (16) n!
Полученное выражение называется первым интерполяционным многочленом Ньютона для интерполирования вперед. Оно справедливо на всем отрезке [x0; xn], но для уменьшения ошибок округления разумно использовать его только для левой половины рассматриваемого отрезка.
Полученная формула для интерполирования вперед практически неудобна для интерполирования функции вблизи конца отрезка. В этом случае используют формулу для многочлена Ньютона для интерполирования назад, которую мы сейчас и получим.
Интерполирующий полином запишем в следующем виде:
N(x) = a0 + a1(x ¡ xn) + a2(x ¡ xn)(x ¡ xn¡1) + :::
+ an(x ¡ xx)(x ¡ xn¡1):::(x ¡ x1): (17) Из условия совпадения значения многочлена и функции в узлах находим:
N(xn) = yn ) a0 = yn;
N(xn¡1) = yn¡1 ) yn + a1(¡h) = yn¡1 )
a1 = yn ¡ yn¡1 = ¢yn¡1 : h h
Аналогично находим
a2 = yn ¡ 2yn¡1 + yn¡2 = ¢2yn¡2 ;
2h2 2h2
и в общем случае, применяя метод математической индукции, можно строго доказать, что
ak = ¢kyn¡k :
k!hk
Подставляя найденные коэффициенты в (17), получаем
N(x) = yn + ¢yn¡1 (x ¡ xn) + 1!h
Делая замену t = x ¡ xn , h
N(x) = N(xn + th) = yn
окончательно находим
+ t¢yn¡1 + t(t + 1)¢2yn¡2 + ::: + 2!
+ t(t + 1):::(t + n ¡ 1)¢ny0: (19) n!
Формула (19) – второй интерполяционный многочлен Ньютона для интерполирования назад.
Отметим, что существует один и только один интерполяционный многочлен при заданном наборе узлов интерполяции. Формулы Лагранжа, Ньютона и др. порождают один и тот же многочлен, разница состоит в алгоритме их построения.
3.3.Точность интерполяции
Значения интерполяционного многочлена y = '(x) и рассматриваемой функции y = f(x) в узлах x = xi, (i = 0; 1; :::; n) в точности совпадают. Если исследуемая функция многочлен степени n, то f(x) ´ '(x). В остальных случаях разность
R(x) = f(x) ¡ '(x) =6 0:
Очевидно, что R(x) есть погрешность интерполяции и называется остаточным членом интерполяционной формулы. Можно показать, что остаточный член интерполяционного многочлена Лагранжа имеет вид
RL(x) = |
(x ¡ x0)(x ¡ x1):::(x ¡ xn) |
f(n+1)(x0): |
(20) |
|
(n + 1)! |
|
В этой формуле f(n+1)(x0) – производная (n + 1)-го порядка функции
f(x) в некоторой точке x0 2 [x0; xn].
Из анализа (20) следует, что погрешность интерполяции тем выше, чем ближе точка x лежит к концам отрезка [x0; xn]. Если же использовать интерполяционный многочлен вне отрезка [x0; xn], то погрешность возрастает очень заметно.
Остаточный член интерполяцинного многочлена Ньютона для случая равноотстоящих узлов следует из (20)
RN (x) = t(t ¡ 1):::(t ¡ n)f(n+1)
(n + 1)!
(x0)hn+1; t = |
x ¡ x0 |
: |
(21) |
|
h |
||||
|
|
|
7
Из вида остаточного члена следует, что повышение степени интерполяционного многочлена уменьшает погрешность, однако из-за неясного поведения f(n+1)(x) возможны проблемы. Поэтому на практике для повышения точности целесообразно уменьшать шаг и выбирать специальное расположение узлов (например, сгущая их к концам отрезка). При этом, как правило, стараются использовать многочлены малой степени.
3.4.Задания для самостоятельной работы
Построить интерполяционный многочлен по заданным таблицам а) в форме Лагранжа; б) в форме Ньютона.
|
xi |
0.2 |
0.5 |
1 |
|
|
1. |
yi |
0.2386693308 |
0.7294255386 |
1.841470985 |
|
|
xi |
1.5 |
2 |
2.2 |
|
|
|
|
|
|
||||
|
yi |
3.247494987 |
4.909297427 |
5.648496404 |
|
|
|
|
|
|
|
|
|
|
xi |
6.0 |
6.3 |
6.6 |
|
|
2. |
yi |
34.286714820 |
39.701203180 |
41.70367307 |
|
|
xi |
6.75 |
6.9 |
7.5 |
|
|
|
|
|
|
||||
|
yi |
41.13764565 |
39.41511178 |
20.43623661 |
|
|
|
|
|
|
|
|
|
|
xi |
2.2 |
2.4 |
2.8 |
|
|
3. |
yi |
2.4114988830 |
2.2626062840 |
2.057777659 |
|
|
xi |
3.2 |
3.4 |
3.8 |
|
|
|
|
|
|
||||
|
yi |
2.001705224 |
2.033201807 |
2.209032288 |
|
|
|
|
|
|
|
|
|
|
xi |
2.0 |
2.2 |
2.4 |
|
|
4. |
yi |
2.9835460860 |
4.0403875470 |
4.756899984 |
|
|
xi |
2.6 |
2.8 |
3.2 |
|
|
|
|
|
|
||||
|
yi |
4.373250602 |
2.640288276 |
-1.283458022 |
|
|
|
|
|
|
|
|
|
|
xi |
0.1 |
0.3 |
0.5 |
|
|
5. |
yi |
9.9833416650 |
3.283557852 |
1.917702154 |
|
|
xi |
0.7 |
1.1 |
1.5 |
|
|
|
|
|
|
||||
|
yi |
1.314729974 |
0.7365350083 |
0.4433311052 |
|
8
|
xi |
4.0 |
4.3 |
4.5 |
|
||
6. |
|
yi |
2.4323981290 |
2.3608144630 |
2.348313255 |
|
|
xi |
4.6 |
4.7 |
4.9 |
|
|||
|
|
||||||
|
|
yi |
2.351940650 |
2.361751112 |
2.398498400 |
|
|
|
|
|
|
|
|
|
|
|
xi |
1.4 |
1.6 |
1.8 |
|
||
7. |
|
yi |
2.4323981290 |
4.2961499960 |
5.844778271 |
|
|
xi |
2.2 |
2.4 |
2.6 |
|
|||
|
|
||||||
|
|
yi |
8.719672866 |
9.428320962 |
9.134725689 |
|
|
|
|
|
|
|
|
|
|
|
xi |
2.0 |
2.3 |
2.9 |
|
||
8. |
|
yi |
3.4864546310 |
4.0469158410 |
5.932133517 |
|
|
xi |
3.2 |
3.8 |
5.0 |
|
|||
|
|
||||||
|
|
yi |
7.354448569 |
11.55532021 |
29.73936426 |
|
|
|
|
|
|
|
|
|
|
|
xi |
0.2 |
0.4 |
0.6 |
|
||
9. |
|
yi |
0.6896558743 |
1.390438928 |
1.774466077 |
|
|
xi |
0.8 |
1.0 |
1.2 |
|
|||
|
|
||||||
|
|
yi |
1.503270954 |
0.3836039536 |
-1.469219613 |
|
|
|
|
|
|
|
|
|
|
|
|
xi |
|
1.2 |
1.25 |
1.35 |
|
10. |
|
yi |
|
4.168650116 |
4.426278450 |
3.896391075 |
|
|
xi |
|
1.40 |
1.45 |
1.50 |
|
|
|
|
|
|
||||
|
|
yi |
|
2.748858752 |
0.7599941296 |
-2.194454593 |
|
3.5.Примеры процедур в среде Maple
3.5.1.Построение интерполяционного многочлена Лагранжа
Для наглядности, построим таблицу значений функции, аналитический вид которой известен.
>restart;
>f:=x-> ln(x)-(sin(x))ˆ 2-exp(x); # аналитический вид функции
>n:=5; # количество шагов
>x0:=1; # первый узел интерполирования
>sh:=0.2; # шаг интерполирования
# задание узлов интерполирования и значений функции в этих узлах > for i from 0 to n do
9
x[i]:=x0+sh*i; y[i]:=f(x[i]); end do;
# Процедура создания интерполяционного многочлена в форме Лагран-
жа
> lagr:=proc(n,x,y,xx)
#n количество узлов интерполяции (степень интерполяционного многочлена)
#x узлы интерполяции
#y значения функции в узлах интерполяции
#xx неизвестная переменная в интерполяционном многочлене local i,j,s,sl;
#i, j переменные циклов
#sl слагаемое в интерполяционном многочлене
#s переменная, для накопления суммы слагаемых (многочлен Лагранжа)
s:=0;
for i from 1 to n do sl:=y[i]; for j from 1 to n do
#получение слагаемых
if (i<>j) then sl:=sl*(xx-x[j])/(x[i]-x[j]); end if; end do;
# накопление всех слагаемых образование многочлена Лагранжа s:=s+sl;
end do; return s; end proc;
> poll:=z-> lagr(n,x,y,z);
# сравнение функции с построенным интерполяционным многочленом
> plot([f(x),poll(x)],x=0..2);
3.5.2.Построение интерполяционного многочлена Ньютона
Для наглядности, построим таблицу значений функции, аналитический вид которой известен.
> restart;
#представление функции в аналитическом виде > yf :=x->sin(x)*ln(x);
#количество шагов, первый узел и шаг интерполирования > N := 20: x0 := 1: h := 0.2:
10