WWW.DISSERS.RU

БЕСПЛАТНАЯ ЭЛЕКТРОННАЯ БИБЛИОТЕКА

   Добро пожаловать!


Pages:     || 2 |
К.П. ЛОВЕЦКИЙ, Л.А. СЕВАСТЬЯНОВ, Е.Б. ЛАНЕЕВ УЧЕБНО-МЕТОДИЧЕСКОЕ ПОСОБИЕ ПО КУРСУ «ВЫЧИСЛИТЕЛЬНЫЙ ЭКСПЕРИМЕНТ И МЕТОДЫ ВЫЧИСЛЕНИЙ» Москва Издательство Российского университета дружбы народов 2007 Ут в е р жд е н о РИС Ученого совета Российского университета дружбы народов Рецензент- кандидат физико-математических наук, доцент Третьяков Н.П.

Ловецкий К.П., Севастьянов Л.А., Ланеев Е.Б.

Учебно-методическое пособие по курсу «Вычислительный эксперимент и методы вычислений». М.: РУДН, 2007. – 35 с.

В пособии изложены примеры программной реализации задач полиномиальной интерполяции. Рассмотрены интерполяционные полиномы Лагранжа, Ньютона и Эрмита. Отдельное внимание уделено исследованию погрешности интерполяции.

Подготовлено на кафедре систем телекоммуникаций.

© Российский университет дружбы народов, Издательство, 2007 ©Ловецкий К.П., Севастьянов Л.А., Ланеев Е.Б., 2007 2 ОГЛАВЛЕНИЕ 1. ИНТЕРПОЛЯЦИОННЫЙ МНОГОЧЛЕН ЛАГРАНЖА 4 Теоретическая часть (реализация примера). 9 Программная реализация примера. 11 2. ИНТЕРПОЛЯЦИОННЫЙ ПОЛИНОМ НЬЮТОНА 16 2.1. Конечные разности 16 2.2.Конечноразностные интерполяционные формулы 19 2.3. Интерполяционная формула Ньютона для неравноотстоящих узлов 20 3. ИНТЕРПОЛЯЦИОННЫЙ ПОЛИНОМ ЭРМИТА 25 3.1. Интерполяция с кратными узлами 25 3.2. Обобщенная интерполяционная формула Ньютона с разделенными разностями 27 4. СХОДИМОСТЬ ИНТЕРПОЛЯЦИОННОГО ПРОЦЕССА 30 5. ИСПОЛЬЗОВАННЫЕ ИСТОЧНИКИ 35 3 Задача аппроксимации функции В основе большинства численных методов математического анализа лежит замена одной функции f (x) (известной, неизвестной или частично известной) другой функцией (x), близкой к f (x) и обладающей «хорошими» свойствами, позволяющими легко производить над нею те или иные аналитические или вычислительные операции.

Будем называть такую замену аппроксимацией.

Будем считать, что аппроксимация функции производится с помощью многочленов степени n N0. Тогда в зависимости от выбора критерия согласия и, в частности, от количества точек согласования f (x) с (x) (будем называть их узлами), то есть точек, в которых известна информация об f (x) и, возможно, ее производных, можно рассмотреть разные конкретные способы аппроксимации.

1. ИНТЕРПОЛЯЦИОННЫЙ МНОГОЧЛЕН ЛАГРАНЖА 1.1. Задача полиномиальной интерполяции Пусть в точках x0, x1,..., xn таких, что a x0 <... < xn b, известны значения функции y = f x, то ( ) есть на отрезке [a,b] задана табличная (сеточная) функция x x0 x1... xn f (x) : (1.1) y y0 y1... yn Функция (x) называется интерполирующей или интерполяционной для f (x) на [a,b], если ее значения x0, x1,..., xn в заданных точках x0, x1,..., xn, назы( ) ( ) ( ) ваемых узлами интерполяции, совпадают с заданными значениями функции f (x), то есть с y0, y1,..., yn соответственно. Тогда задача интерполяции, точнее, полиномиальной, алгебраической или параболической интерполяции (поскольку график любого многочлена называют параболой) формулируется так:

для функции f (x), заданной таблицей (1.1), найти многочлен Pn (x) такой, что выполняется совокупность условий интерполяции Pn(xi ) = yi i 0,1,..., n. (1.2) { } Найти многочлен Pn (x) - это значит, учитывая его каноническую форму Pn (x) = a0 + a1x + a2x2 +...+ anxn, (1.3) найти его n +1 коэффициент a0, a1,...an. Для этого имеется как раз n +1 условие (1.2). Таким образом, чтобы многочлен (1.3) был интерполяционным для функции (1.1) нужно, чтобы его коэффициенты a0, a1,...an удовлетворяли системе уравнений 2 n a0 + a1x0 + a2x0 +...+ anx0 = y 2 n a0 + a1x1 + a2x1 +...+ anx1 = y..............................................

2 n a0 + a1xn + a2xn +...+ anxn = yn.

Из курса алгебры известно, что определитель этой линейной системы (так называемый определитель Вандермонда) отличен от нуля, если все узлы различны, то есть решение этой системы существует и единственно. Однако практическое построение интерполяционного многочлена путем решения такой системы уравнений малоэффективно.

1.2. Базисные полиномы Лагранжа Попробуем построить многочлен n-ой степени Ln(x) n в виде линейной комбинации (x) многочленов одинаcl i i i=ковой n-й степени li (x), i = 0,1,..., n. (Здесь i – номер многочлена от 0 до n.) Для того, чтобы такой многочлен был интерполяционным для функции f x, достаточно зафикси( ) ровать в качестве коэффициентов ci этой линейной комбинации заданные в таблице (1.1) значения yi = f (xi ), а от базисных многочленов li (x) потребовать выполнения условия 0, i j li (xj ) = ij := {} 1, i = j, j,i 0,1,..., n. (1.4) В таком случае для многочлена n Ln (x) = yili (x) i=в каждом узле xj, j 0,1,..., n, в силу (1.4), справедливо { } соотношение Ln(xj ) = l0(xj )y0 +...+ lj-1(xj )yj-1 + lj (xj )yj + lj+1(xj )yj+1 +...

...+ ln(xj )yn = 0 +...+ 0 + yj + 0 +...+ 0 = yj, то есть выполняются условия интерполяции (1.2).

Чтобы конкретизировать вид базисных многочленов li (x), учтем, что они должны удовлетворять условиям (1.4).

Равенство нулю i-го многочлена во всех узлах, кроме i-го, означает, что li (x) можно записать в виде li (x) = Ai (x - x0)...(x - xi-1)(x - xi+1)...(x - xn), a коэффициент Ai этого представления легко получается из содержащегося в (1.4) требования li (xi ) =1. Подставляя в выражение li (x) значение x = xi и приравнивая результат единице, получаем Ai =.

(xi - x0)...(xi - xi-1)(xi - xi+1)...(xi - xn ) Таким образом, базисные многочлены Лагранжа имеют вид (x - x0)(x - x1)...(x - xi-1)(x - xi+1)...(x - xn ) li (x) =, (xi - x0)(xi - x1)...(xi - xi-1)(xi - xi+1)...(xi - xn) а искомый интерполяционный многочлен Лагранжа есть n (x - x0)...(x - xi-1)(x - xi+1)...(x - xn) Ln(x) = yi. (1.5) (xi - x0)...(xi - xi-1)(xi - xi+1)...(xi - xn) i= В качестве примера запишем интерполяционные многочлены Лагранжа первой и второй степени.



При n=1 информация об интерполируемой функции y = f x сосредоточена в двух точках: (x0, y0) и (x1, y1).

( ) Многочлен Лагранжа в этом случае составляется с помощью двух базисных многочленов первой степени (l0(x) и l1(x) ) и имеет вид x - x1 x - xL1(x) = y0 + y1. (1.6) x0 - x1 x1 - xПри n=2 по трехточечной таблице x x0 x1 xf (x) :

y y0 y1 yможно образовать три базисных многочлена (l0(x), l1(x) и l2(x) ), соответственно, интерполяционный многочлен Лагранжа второй степени имеет вид (x - x1)(x - x2) (x - x0)(x - x2) L2(x) = y0 + y1 + (x0 - x1)(x0 - x2) (x1 - x0)(x1 - x2) (1.7) (x - x0)(x - x1) + y2.

(x2 - x0)(x2 - x1) Приближенные равенства f (x) L1(x) и f (x) L2(x) называют соответственно формулами линейной и квадратичной интерполяции.

Пусть для данной функции f x интерполяционный ( ) многочлен Ln(x) построен, то есть для приближенного представления функции f x на отрезке [a,b] [x0, xn] ( ) применяется интерполяционная формула f (x) Ln (x) (1.8) Естественно встает вопрос: какова погрешность такого приближенного равенства Иначе говоря, насколько большим может быть различие между значениями интерполируемой функции f x и соответствующими значениями ( ) интерполяционного многочлена Лагранжа Ln(x) в точках отрезка [a,b], не совпадающих с узловыми точками 1.3. Оценка погрешности интерполяции Знание остаточного члена в предположении (n+1)кратной дифференцируемости f x позволяет записать ( ) точное представление f x через ее интерполяционный ( ) многочлен Ln(x) :

(n+1) f ( ) f (x) = Ln(x) + (x), (1.9) n+(n +1)! где - некоторая (вообще говоря, неизвестная, причем зависящая от x ) точка из промежутка интерполяции (a,b), а n (x) = - xi ) = (x - x0)(x - x1)...(x - xn ) - определен(x n+i=ный через узлы x0, x1,..., xn многочлен (n+1)-ой степени.

Так, если известна величина (n+1) Mn+1 = max f ( ), (1.10) x[a,b] то оценить абсолютную погрешность интерполяционной формулы (1.8) в любой точке x [a,b] можно с помощью неравенства Mn+ Rn(x) = f (x) - Ln (x) n+1(x). (1.11) (n +1)Максимальная погрешность интерполирования на отрезке [a,b] оценивается величиной Mn+ max Rn(x) max n+1(x). (1.12) x[a,b] x[a,b] (n +1)Пример 1.

Теоретическая часть (реализация примера).

Рассмотрим квадратичную интерполяцию функции y = sin x на отрезке 0, / 2 по ее трем значениям:

( ) [ ] sin 0 = 0, sin / 4 = 2 / 2, sin / 2 =1.

( ) ( ) ( ) Для построения полинома Лагранжа нам потребуется три базисных полинома (x - x1)(x - x2) l0(x) =, (x0 - x1)(x0 - x2) (x - x0)(x - x2) l1(x) =, (x1 - x0)(x1 - x2) (x - x0)(x - x1) l2(x) =.

(x2 - x0)(x2 - x1) Полином Лагранжа, согласно формуле (1.7) имеет вид (x - x1)(x - x2) (x - x0)(x - x2) L2(x) = f (x0) + f (x1) + (x0 - x1)(x0 - x2) (x1 - x0)(x1 - x2), (x - x0)(x - x1) + f (x2).

(x2 - x0)(x2 - x1) следовательно, в нашем случае (x - / 4)(x - / 2) (x - 0)(x - / 2) L2(x) =0 + 2 / 2 + (0 - / 4)(0 - / 2) ( / 4 - 0)( / 4 - / 2) (x - 0)(x - / 4) +1, ( / 2 - 0)( / 2 - / 4) или, после преобразований, (1-.

L2(x) = x 2)x + 2 / 2 -1/ 4 (1.13) ( ) Остаточный член для этого случая получаем по формуле (1.11), учитывая, что n = 2, (sin x) ''' =- cos x и (x) = x(x - / 4)(x - / 2).

Имеем:

- cos R2(x) = sin x - L2(x) = x(x - / 4)(x - / 2).

3! Так как точка (0, / 2) неизвестна, можно делать лишь оценки R2(x), полагая M3 = max - cos x = 1. Найдя макx 0, / [ ] симальное значение (x), реализуемое в двух точках данного отрезка x1,2 = (3 + 3) и не превосходящее 0.568, оцениваем сверху величину допустимого отклонения дуги параболы (1.13) от данной синусоиды на промежутке интерполирования:

1 0.max sin x - L2(x) max L2(x) < 0.095.

x 0, / 2 x 0, / [ ] [ ] 3! Подставим в полученный интерполяционный многочлен (1.13) контрольную точку x = / 6. Получим приближенное 4 2 -значение L2( / 6) = 0.517, отличающееся от значения sin( / 6) = 0.5 на величину 0.17, меньшую, чем это допускается оценкой по формуле (1.12):

M3 1 sin( / 6) - L2( / 6) ( / 6) = / 6 0.075.

( ) 3! Наблюдаем типичную картину: фактическая погрешность меньше оценки погрешности в точке, которая, в свою очередь, меньше максимальной погрешности на отрезке (то есть по чебышевской норме).

Программная реализация примера.

Приведем текст подпрограмм-функций на паскале для вычисления значений многочлена Лагранжа.

Подпрограмма-функция вычисления y = sin x :

( ) function f_sinx(x:real):real;

begin result := sin(x);

end;

Процедура вычисления полинома Лагранжа L2(x) (1.13):

function L_2(x:real):real;

begin result := 8/sqr(pi)*x*((1-sqrt(2))*x+(sqrt(2)/2-1/4)*pi);

end;

Процедура оценки погрешности на отрезке 0, / 2 по [ ] nPoints значениям:

function maxDeviation(xLeft, xRight : real;

nPoints : integer) : real;

var i : integer;

x1, s, sw : real;

begin s := 0;

for i:=1 to nPoints do begin x1:=xLeft+(i-1)*(xRight-xLeft)/(nPoints-1);

sw := abs(f_sinx(x1)-L_2(x1));

if s < sw then s := sw;

end;

result := s;

end;

Текст программы и исполняемый модуль находятся на сайте кафедры: http://mathmod.sci.pfu.edu.ru Главное окно с результатом работы программы представлено на рисунке 1. Сплошная линия – график функции y = sin x, штриховая – соответствующий полином Лагранжа (1.13).





Рис1. Интерполяция функции y = sin x Пример 2.

Рассмотрим функцию y = sin 4x на отрезке 0, / 2 и вы[ ] числим полиномы Лагранжа Ln(x).

Приведем текст подпрограмм-функций на паскале для вычисления значений многочлена Лагранжа для произвольного значения n и погрешности интерполяции на заданном количестве точек.

Процедура вычисления полиномов Лагранжа Ln(x) :

type OneArgFunction = function (x:real) : real;

function L_n(x1:real; x : array of real; n:integer; f :

OneArgFunction):real;

var s, sum : real;

i, j : integer;

begin sum:=0;

for i:=1 to n do begin s:=1;

for j:=1 to n do begin if j<>i then s:=s*(x1-X[j])/(X[i]-X[j]) end;

sum:=sum+s*f(X[i]);

end;

result := sum;

end;

Процедура оценки погрешности на отрезке 0, / 2 по [ ] nPoints значениям:

function maxDeviation(xLeft, xRight : real; nPoints : integer) :

real;

var i : integer;

x1, s, sw : real;

begin s := 0;

for i:=1 to nPoints do begin x1:=xLeft+(i-1)*(xRight-xLeft)/(nPoints-1);

sw := abs(f_sinx(x1)-L_2(x1));

if s < sw then s := sw;

end;

result := s;

end;

На рисунке 2 приведен график функции (сплошная линия) и соответствующий ей полином Лагранжа (штриховая линия), построенный по 5-ти точкам. Минимальное отклонение для приведенного метода расчета полинома Лагранжа достигает величины 3.0967E-13 при числе точек интерполяции равном Рис2. Интерполяция функции y = sin 4x по 5-ти точкам Рис3. Интерполяция функции y = sin 4x по 65-ти точкам 20. На рисунке 3 приведен график интерполяционного полинома, построенного по 65-ти точкам на сетке с равномерным шагом. При этом наблюдается значительное расхождение между графиком функции и полиномиальной интерполяцией, особенно заметном на краях интервала.

Образец текста программы и исполняемый модуль находятся на сайте кафедры: http://mathmod.sci.pfu.edu.ru … 2. ИНТЕРПОЛЯЦИОННЫЙ ПОЛИНОМ НЬЮТОНА 2.1. Конечные разности Приведем интерполяционную формулу к более простому виду, подобному виду формулы Тейлора. Если в интерполяционном многочлене Лагранжа (1.5) все слагаемые однотипны и играют одинаковую роль в формировании результата, хотелось бы иметь такое представление интерполяционного многочлена, в котором, как и многочлене Тейлора, слагаемые располагались бы в порядке убывания их значимости. Тогда можно было бы более просто перестраивать его степень, добавляя или отбрасывая удаленные от начала его записи члены.

Рассмотрим сначала частный случай постановки задачи интерполяции. А именно, будем считать, что интерполируемая функция y = f x задана своими значениями ( ) y0, y1,..., yn на системе равноотстоящих узлов x0, x1,..., xn, то есть таких, что любой узел xi этой сетки можно представить в виде xi = x0 + ih, где i = 0,1,..., n, а h > 0 - некоторая постоянная величина, называемая шагом сетки (таблицы).

Прежде чем строить интерполяционные формулы, рассмотрим элементы теории конечных разностей.

Вычитая из каждого последующего члена конечной последовательности из n +1 чисел y0, y1,..., yn предыдущий, образуем n конечных разностей первого порядка.

y0 = y1 - y0, y1 = y2 - y1,...,yn-1 = yn - yn-1, или, проще, n первых разностей данной табличной функции. Из них, в свою очередь, таким же образом можно получить n -1 конечных разностей второго порядка, или вторых разностей:

2 y0 = y1 - y0, 2 y1 = y2 - y1,...,.

2 yn-2 = yn-1 - yn-Этот процесс построения разностей может быть продолжен, и весь он, в результате, описывается одной рекуррентной формулой, выражающей конечную разность k-го порядка k yi через разности (k-1)-го порядка k yi = k -1yi+1 - k -1yi, где k=1,2,…,n и 0 yi = yi.

Имеется прямая связь между конечными разностями и производными. А именно, если учесть, что yi+1 - yi f (xi + h) - f (xi ) ' lim = lim = f (xi ), h0 hhh то можно сказать, что при малых h имеет место приближенное равенство ' yi f (xi )h, то есть первые разности характеризуют первую производную функции f (x), по значениям которой они составлены.

Пользуясь этим, имеем для вторых разностей:

yi+2 - yi+1 yi+1 - yi 2 yi yi+1 - yi h h == h2 h2 h, '' f (xi + h) - f (xi ) '' f (xi ) h '' то есть 2 yi f (xi )h2, и, вообще, (k ) k yi f (xi )hk. (2.1) Таким образом, на конечные разности можно смотреть как на некоторый аналог производных.

Для функции y = f x, заданной таблицей своих зна( ) чений y0, y1,..., yn в узлах x0, x1,..., xn, где xi = x0 + ih, конечные разности разных порядков удобно помещать в одну общую таблицу с узлами и значениями функции (последние можно интерпретировать как конечные разности нулевого порядка). Такую общую таблицу называют таблицей конечных разностей. Заметим, что кроме так называемого диагонального расположения конечных разностей, когда числа в каждом столбце записываются со смещением в полстроки, как это показано в таблице ниже, часто применяют и горизонтальное расположение, где yi, 2 yi и другие разности с индексом i помещают в одной строке с xi, yi.

Таблица 1.

x y 0 yx y 2 y1 y1 3 yx2 y 2 y1 4 y0....

y2 3 yx3 y 2 yyx y 4 2.2.Конечноразностные интерполяционные формулы Пусть функция y = f x задана на сетке равноотстоящих ( ) узлов xi = x0 + ih где i = 0,1,..., n и для нее построена таблица 1 конечных разностей.

Pages:     || 2 |










© 2011 www.dissers.ru - «Бесплатная электронная библиотека»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.