авторефераты диссертаций БЕСПЛАТНАЯ БИБЛИОТЕКА РОССИИ

КОНФЕРЕНЦИИ, КНИГИ, ПОСОБИЯ, НАУЧНЫЕ ИЗДАНИЯ

<< ГЛАВНАЯ
АГРОИНЖЕНЕРИЯ
АСТРОНОМИЯ
БЕЗОПАСНОСТЬ
БИОЛОГИЯ
ЗЕМЛЯ
ИНФОРМАТИКА
ИСКУССТВОВЕДЕНИЕ
ИСТОРИЯ
КУЛЬТУРОЛОГИЯ
МАШИНОСТРОЕНИЕ
МЕДИЦИНА
МЕТАЛЛУРГИЯ
МЕХАНИКА
ПЕДАГОГИКА
ПОЛИТИКА
ПРИБОРОСТРОЕНИЕ
ПРОДОВОЛЬСТВИЕ
ПСИХОЛОГИЯ
РАДИОТЕХНИКА
СЕЛЬСКОЕ ХОЗЯЙСТВО
СОЦИОЛОГИЯ
СТРОИТЕЛЬСТВО
ТЕХНИЧЕСКИЕ НАУКИ
ТРАНСПОРТ
ФАРМАЦЕВТИКА
ФИЗИКА
ФИЗИОЛОГИЯ
ФИЛОЛОГИЯ
ФИЛОСОФИЯ
ХИМИЯ
ЭКОНОМИКА
ЭЛЕКТРОТЕХНИКА
ЭНЕРГЕТИКА
ЮРИСПРУДЕНЦИЯ
ЯЗЫКОЗНАНИЕ
РАЗНОЕ
КОНТАКТЫ


Pages:     | 1 | 2 || 4 | 5 |   ...   | 10 |

«С.П. Шарый Курс ВЫЧИСЛИТЕЛЬНЫХ МЕТОДОВ Курс ВЫЧИСЛИТЕЛЬНЫХ МЕТОДОВ С. П. Шарый Институт вычислительных технологий СО РАН ...»

-- [ Страница 3 ] --

Численным дифференцированием называется процесс нахождения значения производной от функции, использующий значения этой функ ции в некотором наборе точек её области определения. Таким образом, если функция задана таблично, т. е. лишь на конечном множестве зна чений аргумента, либо процедура определения значений этой функции не может быть выписана в виде детерминированной программы, то аль тернатив численному дифференцированию нет. В частности, иногда в виде такого чёрного ящика мы вынуждены представлять вычисление значений функции, аналитическое выражение для которой существует, но является слишком сложным или неудобным для дифференцирова ния первыми двумя способами.

96 2. Численные методы анализа В основе численного дифференцирования лежат различные идеи.

Самая первая состоит в том, чтобы доопределить (восстановить) таб лично заданную функцию до функции непрерывного аргумента, к ко торой уже применима обычная операция дифференцирования. Теория интерполирования, которой посвящены предшествующие параграфы, может оказаться в высшей степени полезной при реализации такого подхода. В частности, таблично заданную функцию можно заменить её интерполяционным полиномом, и его производные считать произ водными рассматриваемой функции. Для этого годится также интер поляция сплайнами или какими-либо другими функциями, а в целом описанный выше подход к численному дифференцированию называют интерполяционным подходом.

2.8а Интерполяционный подход Итак, пусть задан набор узлов x0, x1,..., xn [a, b], т. е. сетка с ша гом hi = xi xi1, i = 1, 2,..., n. Кроме того, заданы значения функции f0, f1,..., fn, такие что fi = f (xi ), i = 0, 1,..., n. Ниже мы рассмотрим простейший вариант интерполяционного подхода, в котором использу ется алгебраическая интерполяция.

Начнём со случая, когда применяется интерполяционный полином первой степени, который мы строим по двум соседним узлам сетки, т. е.

по xi1 и xi, i = 1, 2,..., n:

x xi x xi P1,i (x) = fi1 + fi xi1 xi xi xi fi fi1 fi1 xi fi xi = x+, xi xi1 xi xi где у интерполяционного полинома добавлен дополнительный индекс i, указывающий на ту пару узлов, по которым он построен. Поэтому производная равна fi fi1 fi fi P1,i (x) = =.

xi xi1 hi Это значение можно взять за приближение к производной от рассмат риваемой функции на интервале ]xi1, xi [, i = 1, 2,..., n.

Во внутренних узлах сетки x1, x2,..., xn1, т. е. там, где встре чаются два подинтервала, производную можно брать по любой из воз 2.8. Численное дифференцирование можных формул fi fi f (xi ) fx,i := разделённая разность назад, (2.60) xi xi fi+1 fi f (xi ) fx,i := разделённая разность вперёд. (2.61) xi+1 xi Обе они примерно равнозначны и выбор конкретной из них может быть делом соглашения, удобства или целесообразности. Например, от на правления этой разности может решающим образом зависеть устой чивость разностных схем для численного решения дифференциальных уравнений.

Построим теперь интерполяционные полиномы Лагранжа второй степени по трём соседним точкам сетки xi1, xi, xi+1, i = 1, 2,..., n 1.

Имеем (x xi )(x xi+1 ) (x xi1 )(x xi+1 ) P2,i (x) = fi1 + fi (xi1 xi )(xi1 xi+1 ) (xi xi1 )(xi xi+1 ) (x xi1 )(x xi ) + fi+ (xi+1 xi1 )(xi+1 xi1 ) x2 (xi + xi+1 )x + xi xi+ = fi (xi1 xi )(xi1 xi+1 ) x2 (xi1 + xi+1 )x + xi1 xi+ + fi (xi xi1 )(xi xi+1 ) x2 (xi1 + xi )x + xi1 xi + fi+1.

(xi+1 xi1 )(xi+1 xi ) Поэтому 2x (xi + xi+1 ) 2x (xi1 + xi+1 ) P2,i (x) = fi1 + fi (xi1 xi )(xi1 xi+1 ) (xi xi1 )(xi xi+1 ) 2x (xi1 + xi ) + fi+1.

(xi+1 xi1 )(xi+1 xi ) Воспользуемся теперь тем, что xi xi1 = hi, xi+1 xi = hi+1. Тогда xi+1 xi1 = hi + hi+1, а результат предшествующих выкладок может 98 2. Численные методы анализа быть записан в виде 2x xi xi+ f (x) P2,i (x) = fi hi (hi + hi+1 ) (2.62) 2x xi1 xi+1 2x xi1 xi fi + fi+1.

hi hi+1 hi+1 (hi + hi+1 ) Формула (2.62) может применяться при вычислении значения про изводной в произвольной точке x для случая общей неравномерной сетки. Предположим теперь для простоты, что сетка равномерна, т. е.

hi = h = const, i = 1, 2,..., n. Кроме того, для таблично заданной функции на практике обычно наиболее интересны производные в тех же точках, где задана сама функция, т. е. в узлах x0, x1,..., xn. В точке x = xi из (2.62) получаем для первой производной формулу fi+1 fi f (xi ) f,i = (2.63), x 2h называемую формулой центральной разности. Подставляя в (2.62) ар гумент x = xi1 и сдвигая в получающемся результате индекс на +1, получим 3fi + 4fi+1 fi+ f (xi ).

2h Подставляя в (2.62) аргумент x = xi+1 и сдвигая в получающемся ре зультате индекс на (1), получим fi2 4fi1 + 3fi f (xi ).

2h Займёмся теперь выводом формул для второй производной. Исполь зуя интерполяционный полином второй степени, можно найти:

2 2 f (xi ) P2,i (x) = fi1 fi + fi+1.

hi (hi + hi+1 ) hi hi+1 hi+1 (hi + hi+1 ) В частности, на равномерной сетке с hi = h = const, i = 1, 2,..., n имеем fi1 2fi + fi+ f (xi ) (2.64).

h Естественно, что полученные выражения для второй производной не зависят от аргумента x.

2.8. Численное дифференцирование Несмотря на то, что проведённые выше рассуждения основывались на применении интерполяционного полинома Лагранжа, для взятия производных произвольных порядков на сетке общего вида удобнее ис пользовать интерполяционный полином Ньютона, в котором члены яв ляются полиномами возрастающих степеней.

Выпишем ещё без вывода формулы численного дифференцирова ния на равномерной сетке, полученные по четырём точкам, т. е. с при менением интерполяционного полинома третьей степени: для первой производной f (xi ) (2.65) 11fi + 18fi+1 9fi+2 + 2fi+3, 6h f (xi ) (2.66) 2fi1 3fi + 6fi+1 fi+2, 6h f (xi ) (2.67) fi2 6fi1 + 3fi + 2fi+1, 6h f (xi ) (2.68) 2fi3 + 9fi2 18fi1 + 11fi, 6h для второй производной f (xi ) (2.69) 2fi 5fi+1 + 4fi+2 fi+3, h f (xi ) (2.70) fi1 2fi + fi+1, h f (xi ) (2.71) fi3 + 4fi2 5fi1 + 2fi.

h В формуле (2.70) один из четырёх узлов, по которым строилась фор мула, никак не используется, а сама формула совпадает с формулой (2.64), полученной по трём точкам. Отметим красивую двойственность формул (2.65) и (2.68), (2.66) и (2.67), а также (2.69) и (2.71). Неслуча ен также тот факт, что сумма коэффицентов при значениях функции в узлах во всех формулах равна нулю: он является следствием того, что производная постоянной функции нуль.

В связи с численным дифференцированием и во многих других вопросах вычислительной математики чрезвычайно полезно понятие шаблона формулы, под которым мы будем понимать совокупность охва тываемых этой формулой узлов сетки. Более точно, шаблон формулы 100 2. Численные методы анализа xi2 xi1 xi xi+1 xi+ Рис. 2.11. Шаблон формулы второй разностной производной (2.64).

численного дифференцирования это множество узлов сетки, входя щих в правую часть этой формулы, явным образом либо в качестве ар гументов используемых значений функции. Например, шаблоном фор мулы (2.64) для вычисления второй производной на равномерной сетке fi1 2fi + fi+ f (xi ) h являются три точки xi1, xi, xi+1 (см. Рис. 2.11), в которых должны быть заданы fi1, fi, fi+1. Особенно разнообразны формы шаблонов в случае двух и более независимых переменных.

2.8б Оценка погрешности численного дифференцирования Пусть для численного нахождения k-ой производной функции при меняется формула численного дифференцирования, имеющая шаб лон Ш и использующая значения функции в узлах этого шаблона. Если f (x) дифференцируемая необходимое число раз функция, такая что fi = f (xi ) для всех узлов xi Ш, то какова может быть погрешность вычисления f (k) (x) по формуле ? Вопрос этот можно адресовать как к целому интервалу значений аргумента, так и локально, только к той точке xi, которая служит аргументом левой части формулы численного дифференцирования.

Ранее для погрешности интерполирования функций уже были вы ведены выражения (2.23) и (2.24), и потому заманчивой идеей является их прямое дифференцирование с целью ответа на поставленный выше вопрос. Этот путь оказывается очень непростым, так как применение, к примеру, выражения (2.24) требует достаточной гладкости функции (x), о которой мы можем сказать немногое. Даже если эта гладкость имеется у (x), полученные оценки будут содержать производные (x) и пр., о которых мы знаем ещё меньше. Наконец, шаблон некоторых формул численного дифференцирования содержит меньше точек, чем 2.8. Численное дифференцирование это необходимо для построения интерполяционных полиномов нуж ной степени. Такова, к примеру, формула центральной разности для первой производной или формула для второй производной (2.70), по строенная по четырём точкам на основе полинома 3-й степени. Тем не менее, явные выражения для остаточного члена формул численного дифференцирования на этом пути можно получить методом, который напоминает вывод формулы для погрешности алгебраического интер полирования. Подробности изложены, к примеру, в книгах [18, 56].

Рассмотрим ниже детально более простой и достаточно универсаль ный способ оценивания погрешностей, основанный на разложениях по формуле Тейлора. Суть этого способа заключается, во-первых, в вы писывании по формуле Тейлора разложений для функций, входящих в правую часть формулы численного дифференцирования, и, во-вторых, в аккуратном учёте членов этих разложений с целью получения, по возможности, наиболее точного выражения для ошибки.

Поясним эту методику на примере оценки погрешности для форму лы центральной разности (2.63):

fi+1 fi f (xi ) f,i =.

x 2h Предположим, что f C3 [ xi1, xx+1 ], т. е. функция f трижды непре рывно дифференцируема на интервале между узлами формулы. Под ставляя её в (2.63) и разлагая относительно точки xi по формуле Тейло ра с остаточным членом в форме Лагранжа вплоть до членов второго порядка, получим h2 h3 + f (xi ) + hf (xi ) + f,i = f (xi ) + f ( ) x 2h 2 h2 h f (xi ) hf (xi ) + f (xi ) f ( ) 2 h2 + h = f (xi ) + f ( ) + f ( ), 12 где + и некоторые точки из открытого интервала ] xi1, xi+1 [.

Поэтому h2 + h f,i f (xi ) = f ( ) + f ( ) =, x 12 102 2. Численные методы анализа где = 1 (f ( + ) + f ( ). В целом справедлива оценка M3 f,i f (xi ) h, x в которой M3 = max |f ()| для ] xi1, xi+1 [. То есть, на три жды непрерывно дифференцируемых функциях погрешность вычис ления производной по формуле центральной разности равна O(h2 ) для равномерной сетки шага h.

Определение 2.8.1 Станем говорить, что приближённая формула (численного дифференцирования, интегрирования и т. п.) или же при ближённый численный метод имеют p-ый порядок точности (или по рядок аппроксимации), если на равномерной сетке с шагом h её по грешность является величиной O(hp ), т. е. не превосходит Chp, где C константа, не зависящая от h.

Нередко понятие порядка точности распространяют и на неравно мерные сетки, в которых шаг hi меняется от узла к узлу. Тогда роль величины h играет какой-нибудь характерный размер, описывающий данную сетку, например, h = maxi hi. Порядок точности важная ко личественная мера погрешности формулы или метода, и при прочих равных условиях более предпочтительной является та формула или тот метод, которые имеют более высокий порядок точности. Но сле дует чётко осознавать, что порядок точности имеет асимптотический характер и отражает поведение погрешности при стремлении шагов сетки к нулю. Если этого стремления нет и шаг сетки остаётся доста точно большим, то вполне возможны ситуации, когда метод меньшего порядка точности даёт лучшие результаты, поскольку множитель при hp у него меньше.

Другое необходимое замечение состоит в том, что понятие порядка формулы или метода основывается на сравнении скорости убывания погрешности со скоростью убывания степенных функций 1, x, x2,..., xk,..., то есть существенно завязано на степенную шкалу. Иногда (не слишком часто) эта шкала оказывается не вполне адекватной реально му поведению погрешности.

Пример 2.8.1 Пусть на вещественной оси задана равномерная сетка шага h, включающая в себя узлы 0, ±h, ±2h и т. д. Для функции y = g(x) рассмотрим интерполяцию значения g(0) полусуммой (2.72) g(h) + g(h), 2.8. Численное дифференцирование т. е. простейшим интерполяционным полиномом первой степени по уз лам h и h. Каков будет порядок погрешности такой интерполяции в зависимости от h для различных функций g(x)? 1 1 0 1 1 0 Рис. 2.12. Графики функции y = |x| при 0 1 и 1.

Для функции g(x) = |x|, 0, погрешность интерполяции бу дет, очевидно, равна h, так что её порядок равен. Он может быть нецелым числом (в частности, дробным) и даже сколь угодно малым.

5 Рис. 2.13. График функции y = exp(1/x2 ).

Возьмём в качестве g(x) функцию при x = 0, exp x2, g(x) = при x = 0, 0, известную в математическом анализе как пример бесконечно гладкой, но не аналитической (т. е. не разлагающейся в степенной ряд) функции.

Погрешность интерполяции значения этой функции в нуле с помощью формулы (2.72) равна exp(1/h2 ), при h 0 она убывает быстрее 8 Идея этого примера заимствована из пособия [42], задача 4.2.

104 2. Численные методы анализа любой степени h, так что порядок точности нашей интерполяции ока зывается бесконечно большим. Но такой же бесконечно большой поря док точности интерполирования будет демонстрировать здесь функция y = x2 g(x), хотя для неё погрешность h2 exp(1/h2 ) убывает суще ственно быстрее.

Из выкладок, проведённых для определения погрешности форму лы центральной разности, хорошо видна особенность метода разло жений по формуле Тейлора: его локальный характер, вытекающий из свойств самой формулы Тейлора. Наши построения оказываются при вязанными к определённому узлу (или узлам) сетки, относительно ко торого и следует строить все разложения, чтобы обеспечить взаимные уничтожения их ненужных членов. Как следствие, в этом специальном узле (узлах) мы можем быстро оценить погрешность. Но за пределами этого узла (узлов), в частности, между узлами сетки всё гораздо слож нее и не так красиво, поскольку взаимные уничтожения членов могут уже не происходить.

Какой порядок точности имеют другие формулы численного диф ференцирования?

Методом разложений по формуле Тейлора для дважды гладкой функции f нетрудно получить оценки M2 M |fx,i f (xi )| |fx,i f (xi )| (2.73) h, h, 2 где M2 = max |f ()| по из соответствующего интервала между уз лами. Таким образом, разность вперёд (2.60) и разность назад (2.60) имеют всего лишь первый порядок точности. Отметим, что для дважды непрерывно дифференцируемых функций оценки (2.73) уже не могут быть улучшены и достигаются, к примеру, на функции f (x) = x2.

Конспективно изложим другие результаты о точности формул чис 2.8. Численное дифференцирование ленного дифференцирования:

f (xi ) = 3fi + 4fi+1 fi+2 + O(h2 ), 2h f (xi ) = fi2 4fi1 + 3fi + O(h2 ), 2h f (xi ) = 2fi1 3fi + 6fi+1 fi+2 + O(h3 ), 6h f (xi ) = fi2 6fi1 + 3fi + 2fi+1 + O(h3 ).

6h Оценим теперь погрешность формулы (2.64) для второй производ ной fi1 2fi + fi+ f (xi ) fxx,i =.

h Обозначая для краткости fi = f (xi ) и fi = f (xi ), получим для неё h2 h3 h4 (4) fi hfi + f f ( ) 2fi fxx,i = f+ 2i 6i h2 h2 h3 h4 (4) + + fi + hfi + f+ f+ f ( ) 2i 6i h2 (4) = fi + f ( ) + f (4) ( + ), где, + некоторые точки из открытого интервала ]xi1, xi+1 [, и если f C4 [xi1, xi+1 ], то справедлива оценка M4 |f (xi ) fxx,i | h, где M4 = max |f (4) ()|. Таким образом, порядок точности этой форму лы равен 2 на функциях из C4.

Приведём ещё без вывода результат о погрешности формулы для вычисления второй производной вблизи края сетки (таблицы):

f (xi ) = 2fi 5fi+1 + 4fi+2 fi+3 + O(h2 ), h f (xi ) = fi3 4fi2 + 5fi1 2fi + O(h2 ).

h 106 2. Численные методы анализа Порядок этих формул всего лишь второй, откуда видна роль симмет ричности шаблона в трёхточечной формуле (2.64) с тем же порядком точности.

Что произойдёт, если дифференцируемая функция не будет иметь достаточную гладкость? Тогда мы не сможем выписывать необходимое количество членов разложения по формуле Тейлора, и потому порядок точности формул может быть меньшим, чем для функций с достаточно высокой гладкостью. Это иллюстрирует следующий Пример 2.8.2 Рассмотрим функцию g(x) = x |x|, которую эквива лентным образом можно задать в виде x2, если x 0, g(x) = если x 0.

x, Её график изображён на Рис. 2.14.

Функция g(x) дифференцируема всюду на числовой оси. При x = она имеет производную, равную g (x) = = x |x| + x |x| = |x| + x sgn x = 2 |x|, x |x| а в нуле x |x| g (0) = lim = 0.

x0 x Таким образом, производная g (x) = 2 |x| всюду непрерывна. Но она недифференцируема в нуле, так что вторая производная g (0) уже не существует. Как следствие, g(x) C1, но g(x) C2 на любом интервале, содержащем нуль.

Воспользуемся для численного нахождения производной g (0) фор мулой центральной разности (2.63) на шаблоне с шагом h, симметрич ном относительно нуля:

h2 + h g(h) g(h) h |h| (h) | h| g (0) = = = h.

2h 2h 2h Таким образом, при h 0 приближённое числовое значение производ ной стремится к g (0) = 0 c первым порядком по h, а не вторым, как мы установили это ранее для дважды гладких функций.

2.8. Численное дифференцирование Рис. 2.14. График функции y = x |x|: увидеть разрыв её второй производной в нуле почти невозможно.

2.8в Метод неопределённых коэффициентов это другой подход к по Метод неопределённых коэффициентов лучению формул численного дифференцирования, особенно удобный в многомерном случае, когда построение интерполяционного полинома становится непростым.

Предположим, что задан шаблон из p + 1 штук точек x0, x1,..., xp.

Станем искать приближённое выражение для производной от функции в виде формулы p f (k) (x) (2.74) ci f (xi ), i= которая мотивируется тем обстоятельством, что дифференцирование любого порядка является операцией, линейной по значениям функции.

При этом коэффициенты ci постараемся подобрать так, чтобы эта фор мула являлась точной формулой для какого-то достаточно представи тельного набора функций. Например, в качестве таких пробных функ ций можно взять все полиномы степени не выше заданной либо три гонометричексие полиномы (2.4) какого-то фиксированного порядка и т. п. Рассмотрим ниже подробно случай алгебраических полиномов.

Возьмём f (x) равной последовательным степеням переменной x, т. е. 1, x, x2,..., xq для некоторого фиксированного q. Если форму ла (2.74) обращается в точное равенство на этих пробных функциях, 108 2. Численные методы анализа то с учётом её линейности можно утверждать, что она будет точной для любого алгебраического полинома степени не выше q.

Каждое условие, выписанное для какой-то определённой степени xj, j = 0, 1,..., q, является линейным соотношением для неизвестных ко эффициентов ci, и в целом мы приходим к системе линейных уравнений относительно ci, i = 0, 1,..., p. Для разрешимости этой системы есте ственно взять число неизвестных равным числу уравнений, т. е. p = q.

Получающаяся система линейных уравнений имеет вид c0 + c1 +... + cp = 0, c x + c x +... + c x = 0, 00 11 pp....

..

....

.

....

c0 xk1 + c1 xk1 +... + cp xk1 = 0, p 0 (2.75) c0 xk + c1 xk +... + cp xk = k!, p 0 c0 xk+1 + c1 xk+1 +... + cp xk+1 = (k + 1)! x, p 0.....

..

.....

.

.....

c xp + c xp +... + c xp = p(p 1) · · · (p k + 1) xpk.

00 11 pp В правых частях этой системы стоят k-е производные от 1, x, x2,..., xq, а матрицей системы является матрица Вандермонда вида (2.7), ко торая неособенна для попарно различных узлов x0, x1,..., xp. При этом система линейных уравнений однозначно разрешима относитель но c0, c1,..., cp для любой правой части, но содержательным является лишь случай k p. В противном случае, если k p, правая часть системы (2.75) оказывается нулевой, и, как следствие, система также имеет только бессодержательное нулевое решение. Этот факт имеет интуитивно ясное объяснение: нельзя построить формулу для вычис ления производной k-го порядка от функции, используя значения этой функции не более чем в k точках.

Матрицы Вандермонда в общем случае являются плохообусловлен ными (см. §3.5б). Но на практике решение системы (2.75) вручную или на компьютере обычно не приводит к большим ошибкам, так как порядок системы (2.75), равный порядку производной, бывает, как правило, небольшим. 9 На стр. 84 мы имели возможность обсуждать вопрос о том, каков наивысший порядок производных, всё ещё имеющих содержательный смысл.

2.8. Численное дифференцирование Интересен вопрос о взаимоотношении метода неопределённых коэф фициентов и рассмотренного ранее в §2.8а интерполяционного подхода к численному дифференцированию. К примеру, Ш.Е. Микеладзе в кни ге [53] утверждает, что любая формула численного дифференцирова ния, полученная методом неопределённых коэффициентов, может быть выведена также с помощью интерполяционного подхода, отказывая ме тоду неопределённых коэффициентов в оригинальности. Но нельзя от рицать также, что метод неопределённых коэффициентов конструк тивно проще и технологичнее в применении, и уже только это обстоя тельство оправдывает его существование.

2.8г Полная вычислительная погрешность численного дифференцирования Рассмотрим поведение полной погрешности численного дифферен цирования при расчётах на реальных вычислительных устройствах.

Под полной погрешностью мы понимаем суммарную ошибку числен ного нахождения производной, вызванную как приближённым харак тером самого метода, так и неточностями вычислений на цифровых ЭВМ из-за неизбежных ошибок округления и т. п.

Предположим, к примеру, что первая производная функции вычис ляется по формуле разность вперёд fi+1 fi f (xi ) fx,i =.

h Как мы уже знаем, её погрешность M2 h |fx,i f (xi )|, где M2 = max[a,b] |f ()|. Если значения функции вычисляются с ошибками, то вместо точных fi и fi+1 мы получаем их приближённые значения fi и fi+1, такие что и |fi fi | |fi+1 fi+1 |, где через обозначена предельная абсолютная погрешность вычисле ния значений функции. Тогда в качестве приближённого значения про изводной мы должны взять fi+1 fi f (xi ), h 110 2. Численные методы анализа а предельную полную вычислительную погрешность E(h, ) нахожде ния первой производной функции можно оценить следующим образом:

fi+1 fi f (xi ) E(h, ) = h fi+1 fi fi+1 fi fi+1 fi f (xi ) + h h h (fi+1 fi+1 ) + (fi fi ) M2 h (2.76) + h |fi+1 fi+1 | + |fi fi | M2 h 2 M2 h + = +.

h 2 h Отметим, во-первых, что эта оценка, достижима при подходящем сочетании знаков фигурирующих в неравенствах величин, коль скоро достижимо используемое в преобразованиях неравенство треугольника |a + b| |a| + |b| и достижима оценка погрешности (2.73) для формулы разность вперёд. Во-вторых, оценка не стремится к нулю при умень шении шага h, так как первое слагаемое неограниченно увеличивается при h 0. В целом, функция E(h, ) при фиксированном имеет ми нимум, определяемый условием E(h, ) 2 M2 h 2 M = = + + = 0.

h h h h 2 То есть, оптимальное значение шага численного дифференцирования, при котором достигается минимальная полная погрешность, равно h = 2 (2.77) /M2, и брать меньший шаг численного дифференцирования смысла нет. Зна чение достигающейся при этом полной погрешности есть E(h, ) = 2 M2.

Пример 2.8.3 Пусть в арифметике двойной точности с плавающей точкой, реализованной согласно стандарту IEEE 754/854, численно на ходится производная функции, выражение для которой требует десять 2.8. Численное дифференцирование вычислений, а модуль второй производной ограничен сверху величи ной M2 = 10. Погрешность отдельной арифметической операции мож но считать приближённо равной половине расстояния между соседни ми машинно представимыми числами, т. е. примерно 1016 в районе единицы. Наконец, пусть абсолютная погрешность вычисления функ ции складывается из сумм абсолютных погрешностей каждой опера ции, так что 10 · 1016 = 1015 при аргументах порядка единицы.

Тогда в соответствии с формулой (2.77) имеем h = 2 /M2 = 2 · 108, т. е. брать шаг сетки меньше 108 смысла не имеет.

E(h, ) 0 h Рис. 2.15. Типичный график полной погрешности численного дифференцирования Совершенно аналогичная ситуация имеет место и при использова нии других формул численного дифференцирования. Производная k-го порядка определяется в общем случае формулой вида f (k) (x) = hk (2.78) ci f (xi ) + Rk (f, x), i где ci = O(1). Если эта формула имеет порядок точности p, то её остаточный член оценивается как Rk (f, x) c(x) hp. Этот остаточный член определяет идеальную погрешность численного дифференци рования в отсутствие ошибок вычисления функции, и он неограничен но убывает при h 0.

Но если погрешность вычисления значений функции f (xi ) в узлах 112 2. Численные методы анализа Рис. 2.16. Возмущение функции добавкой sin(nx).

n равна, то в правой части (2.78) возникает ещё член, абсолютная ве личина которого совершенно аналогично (2.76) оценивается сверху как hk |ci |.

i Она неограниченно возрастает при h 0. В целом график полной вы числительной погрешности численного дифференцирования выглядит в этом случае примерно так, как на Рис. 2.15.

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

Отмеченное выше обстоятельство потенциально сколь угодно боль шого возрастания погрешности численного дифференцирования, в дей ствительности, является отражением более глубокого факта некоррект ности задачи дифференцирования, когда её решение не зависит непре рывно от входных данных (см. §1.3). Если f (x) исходная функция, производную которой нам требуется найти, то возмущённая функция f (x) + n sin(nx) при n будет равномерно сходиться к исходной, тогда как её производная f (x) + cos(nx) не сходится к производной f (x) (см. Рис. 2.16). При возмущении исход ной функции слагаемым n sin(n2 x) производная вообще может сколь угодно сильно отличаться от производной исходной функции.

2.9. Алгоритмическое дифференцирование 2.9 Алгоритмическое дифференцирование Пусть u = u(x) и v = v(x) некоторые выражения от переменной x, из которых далее с помощью сложения, вычитания, умножения или деления конструируется более сложное выражение. Напомним правила дифференцирования выражений, образованных с помощью элементар ных арифметических операций:

(u + v) = u + v, (2.79) (2.80) (u v) = u v, (2.81) (uv) = u v + uv, u v uv u (2.82) =.

v v Из них следует, что численное значение производной для сложного вы ражения мы можем найти, зная лишь значения образующих его под выражений и их производных.

Сделанное наблюдение подсказывает идею ввести на множестве пар вида (u, u ), которые составлены из значений выражения и его произ водной, арифметические операции по правилам, следующим из формул (2.79)–(2.82):

(u, u ) + (v, v ) = (u + v, u + v ), (2.83) (2.84) (u, u ) (v, v ) = (u v, u v ), (2.85) (u, u ) · (v, v ) = (uv, u v + uv ), (u, u ) u u v uv (2.86) =,.

(v, v ) v v Первые члены пар преобразуются просто в соответствии с применяе мой арифметической операцией, а операции над вторыми членами пар это в точности копии правил (2.79)–(2.82). Если для заданного выра жения мы начнём вычисления по выписанным формулам (2.83)–(2.86), заменив исходные переменные на пары (x, 1) и константы на пары (c, 0), то на выходе получим пару, состоящую из численного значения выражения и производной от него.

Помимо арифметических операций интересующее нас выражение может содержать вхождения элементарных функций. Для них в соот 114 2. Численные методы анализа ветствии с формулами дифференциального исчисления можем опреде лить действия над парами следующим образом exp (u, u ) = (exp u, u exp u), sin (u, u ) = (sin u, u cos u), (u, u ) = u2, 2uu, (u, u ) = u3, 3u2 u и т.д.

Арифметику пар вида (u, u ) с операциями (2.83)–(2.86) называют дифференциальной арифметикой, а основанный на её использовании способ вычисления значений производных носит название алгоритми ческого дифференцирования. Нередко используют также термин авто матическое дифференцирование.

Строго говоря, мы рассмотрели один из возможных способов ор ганизации алгоритмического дифференцирования, который называют прямым режимом. Существует ещё и обратный режим алгоритмиче ского дифференцирования.

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

Идея алгоритмического дифференцирования может быть распро странена на вычисление разделённых разностей (наклонов) функций, а также на вычисление интервальных расширений производных и на клонов (см., к примеру, [64]).

2.10 Приближение функций 2.10а Обсуждение постановки задачи В этом параграфе мы займёмся задачей приближения функций и подробно рассмотрим случай, когда отклонение функций друг от друга измеряется в так называемой среднеквадратичной норме.

К задаче приближения функций естественно приходят в ситуациях, когда методы интерполирования по различным причинам не удовле 2.10. Приближение функций Рис. 2.17. Интерполяция функции, заданной с погрешностью творяют практику. Эти причины могут носить чисто технический ха рактер. К примеру, гладкость сплайна может оказаться недостаточной, либо его построение слишком сложным. Степень обычного интерпо ляционного полинома может быть неприемлемо высокой для данного набора узлов интерполяции (а высокая степень это трудности при вычислении значений полинома). Но причины отказа от интерполяции могут иметь также принципиальный характер. В частности, это проис ходит, если значения функции в узлах x0, x1,..., xn известны неточно.

В этих условиях целесообразна коррекция самой постановки задачи.

Именно, имеет смысл отказаться от требования, чтобы восстанавли ваемая функция g была точно равна значениям fi в узлах x0, x1,..., xn, допустив, к примеру, для g принадлежности её значений некото рым интервалам, т. е. g(xi ) [f i, f i ], i = 0, 1,..., n, f i f i. Наглядно геометрически это означает построение функции g(x) из заданного класса G, которая в каждом узле сетки xi, i = 0, 1,..., n, проходит через некоторый коридор [f i, f i ], см. Рис. 2.17.

Более общая постановка этой задачи предусматривает наличие неко торой метрики (расстояния), которую мы будем обозначать через dist, и с помощью которой можно измерять отклонение вектора значений (g(x0 ), g(x1 ),..., g(xn )) функции g(x) в узлах сетки от вектора задан ных значений (f0, f1,..., fn ). Фактически, здесь dist это какое-то расстояние на пространстве Rn+1 всех (n + 1)-мерных вещественных векторов. Соответствующая постановка задачи приближения (аппрок 116 2. Численные методы анализа симации) будет звучать тогда так:

Для заданных набора узлов xi, i = 0, 1,..., n, на интервале [a, b] и соответствующих им значений fi, i = 0, 1,..., n, и 0, найти функцию g(x) из класса функций G, такую что dist (f, g), где f = (f0, f1,..., fn ) и g = (g(x0 ), g(x1 ),..., g(xn )).

При этом g(x) называют приближающей (аппроксимирующей) функ цией, Важнейшей модификацией поставленной задачи служит задача наилучшего приближения, когда велчина не фиксируется и ищут при ближающую (аппроксимирующую) функцию g(x), которая доставляет минимум расстоянию dist (f, g).

Согласно классификации §2.1, выписанные выше формулировки яв ляются дискретными вариантами общей задачи о приближении функ ции, в которой дискретный набор узлов x0, x1,..., xn уже не фигу рирует, а отклонение одной функции от другой измеряется на всей области их определения:

Для заданных 0, функции f (x) из F и метрики dist найти функцию g(x) из класса функций G, такую что dist (f, g).

Соответствующая общая формулировка задачи о наилучшем прибли жении ставится так:

Для заданных функции f (x) из класса функций F и метрики dist найти функцию g(x) из класса G, на которой достигается нижняя грань расстояний (2.87) от f (x) до функций из G, т. е. удовлетворяющую условию dist (f, g) = inf hG dist (f, h).

Решение g этой задачи обычно называется наилучшим приближением для f в классе G. Отметим, что в каждом конкретном случае существо вание элемента наилучшего приближения требует отдельного исследо вания.

Хороший и весьма важный в практическом отношении инструмент локального приближения функций формула Тейлора. Но при нело кальном приближении, на широких интервалах, она даёт не очень качественные приближения.

До сих пор ничего не было сказано о выборе классов функций F и G, и в наших формулировках они могут быть весьма произвольными.

Но чаще всего предполагают, что F G, и, кроме того, наделяют F 2.10. Приближение функций и G структурой линейного пространства, определяя на них некоторую норму ·. Именно в ней измеряют отклонение функций (непрерывного или дискретного аргумента) друг от друга, так что dist (f, g) = f g.

Соответственно, в задаче наилучшего приближения функции f ищется такой элемент g G, на котором достигается inf hG f h.

Рассмотренные выше постановки задач дают начало большим и важным разделам математики, в совокупности образующим теорию приближения функций (или теорию аппроксимации). Её ветвью явля ется, в частности, теория равномерного приближения, когда отклоне ние функций оценивается в норме f = maxx[a,b] |f (x)| (см. [46, 58]).

Выбор различных норм (т. е. различных мер отклонения функций друг от друга) и различных классов функций обуславливает огромное раз нообразие задач теории приближения.

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

Предложение 2.10.1 Пусть X нормированное линейное простран ство, а U его конечномерное линейное подпространство. Тогда для для любого f X существует элемент наилучшего приближения u U.

Элемент наилучшего приближения, вообще говоря, может быть не единственным. Но при определённых условиях мы можем гарантиро вать его единственность, опираясь лишь на свойства пространства X.

Нормированное пространство X с нормой · называют строго нормированным, если для произвольных x, y X из равенства x+y = x + y следует существование такого скаляра R, что y = x.

Пример 2.10.1 Строго нормированным пространством является R с евклидовой нормой x 2 = x2 + x2. Но нормы x 1 = |x1 | + |x2 | 1 и x = max{|x1 |, |x2 |} на R2, которые эквивалентны норме · (см. §3.3б), не порождают строго нормированное пространство.

Предложение 2.10.2 Пусть X строго нормированное линейное пространство, а U его конечномерное линейное подпространство.

118 2. Численные методы анализа Тогда для для любого f X элемент наилучшего приближения u U единствен.

Доказательство. Предположим, что для элемента f существуют два наилучших приближения m m u = u i u = u i, и (2.88) i i i=1 i= которые определяются наборами коэффициентов (u, u,..., u ) и (u, 1 2 m u,..., u ) разложения по базису i, i = 1, 2,..., m. Тогда 2 m m m u i u i f f = µ 0.

= i i i=1 i= Возьмём середину отрезка прямой, соединяющей u и u, т. е. точку, у которой компоненты разложения по векторам i, равны 1 (u + u ).

i i Имеем m m m 1 u ) i u u i f f f 2 (ui + = i + i i i 2 i=1 i=1 i= m m 1 u u i = µ.

f f i + i i 2 i=1 i= Строгого неравенства здесь быть не может, что очевидно для µ = 0, а для µ 0 означало бы существование элемента, приближающего f лучше, чем два элемента наилучшего приближения u и u. Поэтому необходимо должно выполняться равенство m m 1 u i u i f f + i i 2 i=1 i= m m 1 u i + u i.

f f = i i 2 i=1 i= Но если рассматриваемое пространство строго нормированное, то m m u i = f u i (2.89) f i i i=1 i= 2.10. Приближение функций для некоторого вещественного.

В случае, когда = 1, получаем m (u u ) i, · f= i i 1 i= т. е. f точно представляется в виде линейной комбинации базисных век торов i. Тогда µ = 0, и в силу единственности разложения по базису u = u, откуда f = 0, что соответствует тривиальному случаю.

i i В остающемся случае = 1 для выполнения равенства (2.89) необ ходимо u = u, и тогда два элемента наилучшего приближения (2.88) i i также должны совпадать.

2.10б Задача приближения в евклидовом пространстве Рассмотрим подробно важный частный случай задачи о наилучшем приближении (2.87), в котором • класс F линейное нормированное пространство функций, на котором задано скалярное произведение ·, ·, и с его помощью норма в F определяется как f = f, f, • класс функций G F, из которого выбирается искомый эле мент наилучшего приближения, является конечномерным под пространством в F.

Напомним, что конечномерные линейные векторные пространства, в которых определено скалярное произведение, называются евклидо выми пространствами. Бесконечномерные линейные векторные про странства со скалярным произведением называются гильбертовыми пространствами при дополнительном условии полноты, т. е. существо вания в них предела всякой фундаментальной последовательности от носительно нормы, порождённой этим скалярным произведением. Гиль бертовы пространства являются ближайшим обобщением пространств с привычной нам геометрией.

В условиях постановки задачи, описанной в начале раздела, будем предполагать, что известен {j }m базис m-мерного линейного под j= пространства G F. Мы ищем приближение g для элемента f F в 120 2. Численные методы анализа виде m (2.90) g= cj j.

j= где cj, j = 1, 2,..., m неизвестные коэффициенты, подлежащие опре делению. Если через обозначить квадрат нормы отклонения f от g, то имеем = f g = f g, f g = f, f 2 f, g + g, g m m m (2.91) = f, f 2 cj f, j + cj ck j, k.

j=1 j=1 k= Как видим, есть квадратичная форма от аргументов c1, c2,..., cm плюс ещё некоторые линейные члены относительно cj и постоянное слагаемое f, f. Так как для всех значений аргументов функция принимает только неотрицательные значения, то ясно, что она должна достигать своего минимума.

Действительно, после приведения к главным осям квадратичная форма в составе обязательно должна получить вид суммы квадра тов с положительными коэффициентами, так как иначе вся была бы неограниченной снизу. Но сумма квадратов неограниченно возрастает при увеличении расстояния от нуля до аргумента c = (c0, c1,..., cm ) и растёт быстрее линейных членов. Следовательно, при достаточно боль ших значениях c2 +... + c2 значение также сколь угодно велико.

m Более точно, для любого K существует такое C, что (x) K при c2 +... + c2 C. Получается, что минимум функции не может m достигаться на бесконечности. Минимум непрерывной функции находится где-то в компактном шаре c2 +... + c2 C, т. е. действи m тельно достигается этой функцией.

Для определения минимума функции продифференцируем её по cj, j = 1, 2,..., m, и приравняем полученные производные к нулю:

m (2.92) = 2 f, j + 2 ck j, k = 0.

cj k= Множитель 2 при сумме всех ck j, k появляется оттого, что в двой ной сумме из выражения (2.91) слагаемое с cj возникает дважды, один раз с коэффициентом j, k, а другой раз с коэффициентом k, j.

2.10. Приближение функций В целом, из равенств (2.92) для определения cj получаем систему линейных алгебраических уравнений m (2.93) j, k ck = f, j, j = 1, 2,..., m.

k= Матрица её коэффициентов 1, 1 1, 2... 1, m, 2, 2... 2, n 21 (2.94) (1, 2,..., m ) =,...

..

...

.

...

m, 1 m, 2... m, m называется, как известно, матрицей Грама системы векторов 1, 2,..., m. Из курса линейной алгебры и аналитической геометрии чита телю должно быть известно, что матрица Грама это симметричная матрица, неособенная тогда и только тогда, когда векторы 1, 2,..., m линейно независимы (см., к примеру, [33]). При выполнении этого условия матрица Грама является ещё и положительно определённой.

Таким образом, решение задачи наилучшего среднеквадратичного при ближения существует и единственно в случае, когда 1, 2,..., m образуют базис в подпространстве G.

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

Наиболее простой вид матрица Грама имеет в случае, когда базис ные функции j ортогональны друг другу, т. е. когда j, k = 0 при j = k. При этом система линейных уравнений (2.93) становится диаго нальной и решается тривиально f, j (2.95) cj =, j = 1, 2,..., m.

j, j Соответствующее наилучшее приближение имеет вид суммы m g= cj j j= 122 2. Численные методы анализа и, как известно, называется (конечным) рядом Фурье для f по орто гональной системе векторов {j }m. Коэффициенты (2.95) называют j= при этом коэффициентами Фурье разложения функции f.

Кроме того, в случае ортогонального и близкого к ортогонально му базиса {j }m решение системы (2.93) устойчиво к возмущениям j= в правой части и неизбежным погрешностям вычислений. Но в общем случае базис линейного подпространства G может сильно отличаться от ортогонального, и тогда свойства системы уравнений (2.93) могут быть плохими в том смысле, что её решение будет чрезвычайно чув ствительным к возмущениям и погрешностям.

2.10в Среднеквадратичное приближение функций В этом разделе мы применим развитый выше общий подход к кон кретной задаче наилучшего среднеквадратичного приближения функ ций, заданных на интервале вещественной оси.

Приближение функций в норме, порождённой скалярным произве дением, часто называют среднеквадратичным приближением или про сто квадратичным. Дело в том, что в конечномерной ситуации, когда имеем f = (f0, f1,..., fn ) и g = (g0, g1,..., gn ), обычно n (2.96) f, g = i f i g i, n+1 i= для некоторого положительного весового вектора = (0, 1,..., n ), i 0. То есть, соответствующая норма · такова, что расстояние одной функции до другой есть 1/ n 1 (2.97) f g i (fi gi ) dist (f, g) = =, n+1 i= усреднение квадратов разностей компонент с какими-то весовыми множителями i, i = 0, 1,..., n. В частности, если известна информа ция о точности задания отдельных значений функции fi, то веса i можно назначать так, чтобы отразить величину этой точности, сопо ставляя значениям fi с бльшей точностью бльший вес. Нормирую о о щий множитель n+1 при сумме в (2.96) удобно брать для того, чтобы 2.10. Приближение функций с ростом размерности n (при росте количества наблюдений, измельче нии сетки и т. п.) ограничить рост величины скалярного произведения и обеспечить соизмеримость результатов при различных n.

Если f и g функции непрерывного аргумента, то обычно полагают скалярное произведение равным b (2.98) f, g = (x)f (x)g(x) dx, a для некоторой весовой функции (x) 0. Это выражение с точностью до множителя можно рассматривать как предел (2.96) при n, так как в (2.96) легко угадываются интегральные суммы Римана для интеграла (2.100) по интервалу [a, b] единичной длины при его равно мерном разбиении на подинтервалы. Тогда аналогом (2.97) является расстояние между функциями 1/ b (2.99) f g = (x) f (x) g(x) dist (f, g) = dx.

a В связи с решением рассматриваемой задачи приближения функций часто используют термин метод наименьших квадратов. Фактически, это общее название целого семейства идейно близких методов постро ения приближений, которые основаны на минимизации суммы квад ратов отклонений компонент исходного вектора от приближающего.

В случае, когда рассматривается приближение функций (или вообще элементов) каких-то абстрактных пространств, естественным обобще нием минимизации суммы квадратов является нахождение минимума нормы, порождённой скалярным произведением.

Пример 2.10.2 В качестве примера практического возникновения задачи среднеквадратичного приближения рассмотрим тепловое дей ствие тока I(t) в проводнике сопротивлением R. Мгновенная тепло вая мощность, как известно из теории электричества, равна при этом I 2 (t)R, а полное количество теплоты, выделившееся между моментами времени a и b, равно b I 2 (t)R dt.

a Если мы хотим, скажем, минимизировать тепловыделение рассматри ваемого участка электрической цепи, то нам нужно искать такой ре 124 2. Численные методы анализа Рис. 2.18. Различие равномерного и интегрального (в частности, среднеквадратичного) отклонений функций друг от друга жим её работы, при котором достигался бы минимум выписанного ин теграла, т. е. среднеквадратичное значение тока. В электротехнике его называют действующим или эффективным значением силы тока.

Линейным пространством F, элементы которого будут приближать ся, выступит пространство всех функций, квадрат (т. е. степень 2) ко торых интегрируем на заданном интервале [a, b]. Его называют про странством L2 [a, b], и нам сначала требуется показать, что оно в самом деле является линейным.

Ясно, что если f L2 [a, b], то для любого скаляра c функция cf (x) также интегрируема с квадратом на [a, b]. Далее, с силу очевидного неравенства 2 2 |f (x) g(x)| f (x) + g(x) из интегрируемости второй степени функций f (x) и g(x) следует инте грируемость их произведения на [a, b]. Поэтому существует интеграл b b b b 2 2 f (x) + g(x) dx = f (x) dx + 2 f (x)g(x) dx + g(x) dx, a a a a т. е. сумма f (x) + g(x) также имеет интегрируемый квадрат. Это за вершает доказательство линейности пространства L2 [a, b]. Скалярное произведение в нём задаётся выражением b (2.100) f, g = f (x)g(x) dx, a аналогичным (2.98). В курсах функционального анализа показывает ся, что если интегрирование понимается в смысле Лебега, то L2 [a, b] 2.10. Приближение функций гильбертово пространство, т. е. дополнительно обладает свойством полноты. По этой причине оно очень популярно в самых различных математических дисциплинах, от теории уравнений в частных произ водных до статистики.

Пример 2.10.3 Рассмотрим задачу о среднеквадратичном приближе нии функций из L2 [a, b] полиномами фиксированной степени m. Ска лярное произведение определяется посредством (2.100), а нормой берём f = f (x) dx.

Соответственно, расстояние между функциями определяется тогда как 1/ f g f (x) g(x) dist (f, g) = = dx.

Если в качестве базиса в линейном подпространстве полиномов мы возьмём последовательные степени x2, xm, 1, x,..., то на месте (i, j) в матрице Грама (2.94) размера (m+1)(m+1) будет стоять элемент xi+j1 i1 j x x dx = =, i, j = 1, 2,..., m + i+j1 i+j 0 (сдвиг показателей степени на (1) вызван тем, что строки и столбцы матрицы нумеруются, начиная с единицы, а не с нуля, как последова тельность степеней). Матрица H = (hij ) с элементами hij = 1/(i+j 1), имеющая вид 1 1...

2 3 m 1 1... m+ 2 3 1 1 1... m+2, 3 4....

..

....

.

....

1 1 1... 2m+ m m+1 m+ называется матрицей Гильберта, и она является исключительно плохо обусловленной матрицей (см. §3.5б). Иными словами, решение СЛАУ с этой матрицей является очень непростой задачей.

126 2. Численные методы анализа Пример 2.10.4 Пусть k и l натуральные числа. Поскольку sin(kx) cos(lx) dx = 0, для любых k, l, и 2 sin(kx) sin(lx) dx = 0, cos(kx) cos(lx) dx = 0, 0 для k = l, то базис из тригонометрических полиномов вида 1, cos(2kx), sin(2kx), k = 1, 2,..., является ортогональным на [0, 1] относительно скалярного произведе ния (2.100) с весом (x) = 1. Иными словами, этот базис очень хорош в вычислительном отношении для построения среднеквадратичных при ближений.

Более детальный теоретический анализ и практический опыт пока зывают, что в методе наименьших квадратов в качестве базиса 1, 2,..., n линейного подпространства G F имеет смысл брать системы элементов, ортогональных, возможно, по отношению к другому ска лярному произведению, так как это служит гарантией разумной ма лости внедиагональных элементов матрицы Грама и, как следствие, её не слишком плохой обусловлености.

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

Среднеквадратичные приближения и метод наименьших квадра тов для решения переопределённых систем линейных алгебраических уравнений, которые возникают в связи с задачами обработки наблюде ний, был почти одновременно предложен на рубеже XVIII–XIX веков К.Ф. Гауссом и А.М. Лежандром, причём последний дал новому подхо ду современное название. На практике метод наименьших квадратов очень часто применяется в силу двух главных причин. Во-первых, его применение бывает вызвано ясным содержательным смыслом задачи, в которой в качестве меры отклонения возникает именно сумма квад ратов или интеграл от квадрата функции. К примеру, чрезвычайно 2.11. Полиномы Лежандра популярно теоретико-вероятностное обоснование метода наименьших квадратов, данное впервые также К.Ф. Гауссом и доведённое до со временного состояния в трудах П.С. Лапласа, П.Л. Чебышёва и потом А.А. Маркова и А.Н. Колмогорова. Во-вторых, в методе наименьших квадратов построение элемента наилучшего приближения сводится к решению системы линейных уравнений, т. е. хорошо разработанной за даче. Если для измерения расстояния между функциями применяются какие-то другие метрики, отличные от (2.99), то решение задачи мини мизации этого расстояния может быть существенно более трудным. В целом, если какое-либо одно или оба из выписанных условий не выпол няется, то метод наименьших квадратов может быть не самой лучшей возможностью решения задачи приближения.


Нередко форма приближающей функции (2.90) не подходит по тем или иным причинам, и тогда приходится прибегать к нелинейному ме тоду наименьших квадратов, когда приближающая функция g(x) вы ражается нелинейными образом через базисные функции 1 (x), 2 (x),..., m (x). Тогда минимизация средневадратичного отклонения f от g уже не сводится к решению системы линейных алгебраических урав нений (2.93), и нам нужно применять численные методы оптимизации.

Обсуждение этого круга вопросов и дальнейшие ссылки можно найти, к примеру, в книге [41].

2.11 Полиномы Лежандра 2.11а Мотивация и определение Примеры 2.10.2 и 2.10.3 из предшествующего раздела показывают, что выбор хорошего, т. е. ортогонального или почти ортогонального, ба зиса для среднеквадратичного приближения функций является нетри виальной проблемой. Для её конструктивного решения можно восполь зоваться, к примеру, известным из курса линейной алгебры процессом ортогонализации Грама-Шмидта или его модификациями (см. §3.7е).

Напомним, что по данной конечной линейно независимой системе век торов v1, v2,..., vn этот процесс строит ортогональный базис q1, q2,..., qn линейной облочки векторов v1, v2,..., vn. Он имеет следующие 128 2. Численные методы анализа расчётные формулы:

(2.101) q1 v1, k vk, qi (2.102) qk vk qi, k = 2,..., n.

qi, qi i= Иногда получающийся ортогональный базис дополнительно нормиру ют.

Конкретный вид ортогональных алгебраических полиномов, кото рые получатся в этом процессе из канонического базиса 1, x, x2,..., зависит, во-первых, от интервала [a, b], для которого рассматривается скалярное произведение (2.100), и, во-вторых, от весовой функции (x).

Для частного случая единичного веса, когда (x) = 1, мы можем су щественно облегчить свою задачу, если найдём семейство ортогональ ных полиномов для какого-нибудь одного, канонического, интервала [, ], а затем для любого другого интервала будем пользоваться фор мулой линейной замены переменной y = kx+ l со специальным образом подобранными константами k, l R, k = 0. Тогда x = (y l)/k, и для a = k + l, b = k + l имеем равенство b yl yl f (x) g(x) dx = f g dy, k k k a справедливое в силу формулы замены переменных в определённом ин теграле. Из него вытекает, что равный нулю интеграл по канониче скому интервалу [, ] останется нулевым и при линейной замене пе ременных. Как следствие, получающиеся при такой замене полиномы 1 f ( k (y l)) и g( k (y l)) будут ортогональны на [a, b].

В качестве канонического интервала [, ] обычно берётся [1, 1], и тогда формула замены переменных принимает вид y = 1 (b a) x + 1 (a + b), 2 так что переменная y пробегает интервал [a, b], если x [1, 1]. Обрат ное преобразование даётся формулой 2y (a + b), x= ba которая позволяет построить ортогональные полиномы для любого ин тервала, зная их для [1, 1]. Рассмотренный рецепт может быть при менён не только к построению ортогональных полиномов, но также к любым произвольным функциям.

2.11. Полиномы Лежандра Полиномами Лежандра называют семейство полиномов Ln (x), за висящих от неотрицательного целого параметра n, которые образуют ортогональную систему относительно скалярного произведения (2.100) с простейшим весом (x) = 1 на интервале [1, 1]. Они были введены в широкий оборот французским математиком А. Лежандром в 1785 году.

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

Применяя к степеням 1, x, x2, x3,... последовательно формулы ортогонализации (2.101)–(2.102) со скалярным произведением (2.100) на интервале [1, 1], получим x2 1, x3 (2.103) 1, x, x,...

3 (два первых полинома оказываются изначально ортогональными).

Часто в качестве альтернативного представления для полиномов Лежандра используют формулу Родрига dn 1 n (2.104) x Ln (x) =, n = 0, 1, 2,....

2n n! dxn Очевидно, что функция Ln (x), определяемая этой формулой, является алгебраическим полиномом n-ой степени со старшим коэффициентом, не равным нулю, так как при n-кратном дифференцировании полинома (x2 1)n = x2n nx2(n1) +... + (1)n степень понижается в точности на n. Коэффициент 1/(2n n!) перед производной в (2.104) взят с той целью, чтобы удовлетворить условию Ln (1) = 1. Всюду далее посред ством Ln (x) мы будем обозначать полиномы Лежандра, определяемые формулой (2.104).

Предложение 2.11.1 Полиномы Ln (x), n = 0, 1,..., задаваемые фор мулой Родрига (2.104), ортогональны друг другу в смысле скалярного произведения на L2 [1, 1] с единичным весом. Более точно, 0, если m = n, 1 Lm (x)Ln (x) dx =, если m = n.

2n + Доказательство. Обозначая (x) = (x2 1)n, 130 2. Численные методы анализа можно заметить, что dk 2 n (k) (x) = при x = ±1, x 1 =0 k = 0, 1, 2,..., n 1.

dxk Кроме того, в силу формулы Родрига (2.104) (n) (x), Ln (x) = n = 0, 1, 2,....

2n n!

Поэтому, если Q(x) является n раз непрерывно дифференцируемой функцией на [1, 1], то, последовательно применяя формулу интегри рования по частям, получим 1 Q(x) (n) (x) dx Q(x) Ln (x) dx = 2n n!

1 1 1 Q(x) (n1) (x) Q (x) (n1) (x) dx = n n! n n!

2 2 Q (x) (n1) (x) dx = 2n n! ··· = = (1)n Q(n) (x) (x) dx. (2.105) 2n n! Если Q(x) любой полином степени меньше n, то его n-ая производная Q(n) (x) равна тождественному нулю, а потому из полученной формулы тогда следует Q(x) Ln (x) dx = 0.

В частности, это верно и в случае, когда вместо Q(x) берётся поли ном Lk (x) степени k, меньшей n, что доказывает ортогональность этих полиномов с разными номерами.

Найдём теперь скалярное произведение полинома Лежандра с са мим собой. Если Q(x) = Ln (x), то d2n 1 (2n)!

n Q(n) (x) = x 1 =.

2n n! dx2n 2n n!

2.11. Полиномы Лежандра По этой причине из (2.105) следует 1 (2n)!

Ln (x) Ln (x) dx = (1)n (x) dx (2n n!) 1 (2n)!

(1 x2 )n dx.

= (2n n!)2 Если обозначить (1 x2 )n dx, Zn = то нетрудно найти 1 (1 x2 )n1 dx x2 (1 x2 )n1 dx Zn = 1 x2 (1 x2 )n1 dx.

= Zn Интегрируя по частям вычитаемое в последнем выражении, получим 1 x2 (1 x2 )n1 dx = x d(1 x2 )n 2n 1 1 1 1 x(1 x2 )n (1 x2 )n dx = = + Zn.

2n 2n 2n Комбинируя эти результаты, будем иметь Zn = Zn1 2n Zn, откуда 2n Zn = Zn1.

2n + Для нахождения числовых значений Zn учтём, что 1 (1 x2 )0 dx = Z0 = 1 dx = 2.

1 Тогда Z1 = 2 · и т. д., так что 2 · 4 ·... · 2n Zn = 2.

3 · 5 ·... · (2n + 1) 132 2. Численные методы анализа Окончательно, скалярное произведение Ln (x) на себя равно 1 · 2 ·... · 2n (2n)! Zn = Zn =.

(2n n!)2 (2 · 4 ·... · 2n)2 2n + Это завершает доказательство предложения.

2.11б Основные свойства полиномов Лежандра Выпишем первые полиномы Лежандра, как они даются формулой Родрига (2.104):

L0 (x) = 1, L1 (x) = x, L2 (x) = 2 (3x2 1), L3 (x) = 1 (5x3 3x), (2.106) L4 (x) = 1 (35x4 30x2 + 3), L5 (x) = 8 (63x5 70x3 + 15x), ···.

Они с точностью до множителя совпадают с результатом ортогонали зации Грама-Шмидта (2.103). Графики полиномов (2.106) изображены на Рис. 2.19, и они похожи на графики полиномов Чебышёва. В од ном существенном моменте полиномы Лежандра всё же отличаются от полиномов Чебышёва: абсолютные значения локальных минимумов и максимумов на [1, 1] у полиномов Лежандра различны и не могут быть сделаны одинаковыми ни при каком масштабировании.

Тем не менее, сходство полиномов Лежандра и полиномов Чебышё ва подтверждает следующее Предложение 2.11.2 Все нули полиномов Лежандра Ln (x) вещест венны, различны и находятся на интервале [1, 1].

Доказательство. Предположим, что среди корней полинома Ln (x), лежащих на [1, 1], имеется s штук различных корней 1, 2,..., s 2.11. Полиномы Лежандра 0. 0. 0. 0. 0. 0. 0. 0. 1 0.8 0.6 0.4 0.2 0 0.2 0.4 0.6 0.8 Рис. 2.19. Графики первых полиномов Лежандра на интервале [1.2, 1.2].

нечётной кратности 1, 2,..., s, так что Ln (x) = (x 1 )1 (x 2 )2 · · · (x s )s (x), где в полиноме (x) присутствуют корни Ln (x), не лежащие на [1, 1], а также те корни Ln (x) из [1, 1], которые имеют чётную кратность.

Таким образом, (x) уже не меняет знака на интервале [1, 1]. Ясно, что s n, и наша задача установить равенство s = n.

Рассмотрим интеграл I= Ln (x) (x 1 )(x 2 ) · · · (x s ) dx (x 1 )1 +1 (x 2 )2 +1 · · · (x s )s +1 (x) dx.

= Теперь 1 +1, 2 +1,..., s +1 чётные числа, так что подинтегральное выражение не меняет знак на [1, 1]. Это выражение равно нулю лишь в конечном множестве точек, и потому определённо I = 0.

С другой стороны, выражение для I есть скалярное произведение, в смысле L2 [1, 1], полинома Ln (x) на полином (x1 )(x2 ) · · · (xs ) 134 2. Численные методы анализа степени не более n 1, если выполнено условие s n. Следовательно, в силу свойств полиномов Лежандра при этом должно быть I = 0.

Полученное противоречие может быть снято только в случае s = n, т. е. когда все корни полинома Ln (x) различны и лежат на интервале [1, 1].

Отметим, что проведённое доказательство проходит для скалярных произведений вида (2.98) с достаточно произвольными весовыми функ циями (x), а не только для единичного веса. Кроме того, тот факт, что интервал интегрирования есть [1, 1], также нигде не использовался в явном виде. Фактически, это доказательство годится даже для беско нечных пределов интегрирования. Оно показывает, что корни любых ортогональных полиномов вещественны и различны.

Можно показать дополнительно, что нули полинома Лежандра Ln (x) перемежаются с нулями полинома Ln+1 (x). Кроме того, аналогично по линомам Чебышёва, нули полиномов Лежандра сгущаются к концам интервала [1, 1].


Ещё одно интересное свойство полиномов Лежандра, задаваемых посредством формулы Родрига (2.104):

Ln (1) = (1)n, Ln (1) = 1, n = 0, 1, 2,....

Кроме того, справедливо рекуррентное представление (n + 1)Ln+1 (x) = (2n + 1) xLn (x) nLn1 (x).

Доказательство этих свойств можно найти, в частности, в [26, 53]. По следняя формула даёт практически удобный способ вычисления зна чений полиномов Лежандра, так как в их явном представлении (2.106) коэффициенты растут экспоненциально быстро в зависимости от номе ра полинома и, как следствие, прямые вычисления с ними могут дать большую погрешность.

Введём так называемые приведённые полиномы Лежандра Ln (x), старший коэффициент у которых равен единице. Как следствие фор мулы Родрига (2.104) можем выписать следующее представление:

dn 2 dn 1 n!

n n x 1 x Ln (x) = =, 2n (2n 1) · · · (n + 1) dxn (2n)! dxn n = 1, 2,.... Как и исходная формула Родрига, выражение после вто рого равенства имеет также смысл при n = 0, если под производной нулевого порядка от функции понимать её саму.

2.11. Полиномы Лежандра Предложение 2.11.3 Среди всех полиномов степени n, n 1, со старшим коэффициентом, равным 1, полином Ln (x) имеет на ин тервале [1, 1] наименьшее среднеквадратичное отклонение от нуля.

Иными словами, если Qn (x) полином степени n со старшим коэф фициентом 1, то 1 2 (2.107) dx Qn (x) Ln (x) dx.

1 Доказательство. Если Qn (x) = xn + an1 xn1 +... + a1 x + a0, то для отыскания наименьшего значения выражения J (a0, a1,..., an1 ) = Qn (x) dx (2.108) (xn + an1 xn1 +... + a1 x + a0 )2 dx = продифференцируем этот интеграл по коэффициентам a0, a1,..., an и приравняем полученные производные к нулю. Так как в данных усло виях дифференцирование интеграла по параметру, от которого зависит подинтегральная функция, сводится к взятию интеграла от производ ной, то имеем в результате J 2 xn + an1 xn1 +... + a1 x + a0 xk dx = ak Qn (x) xk dx = 0, = k = 0, 1,..., n 1. Это означает, что полином Qn (x) ортогонален в смысле L2 [1, 1] всем полиномам меньшей степени. Следовательно, при минимальном значении интеграла (2.108) полином Qn (x) обязан совпа дать с полиномом Лежандра.

Для построения полинома, который имеет наименьшее среднеквад ратичное отклонение от нуля на произвольном интервале [a, b] можно воспользоваться линейной заменой переменной с масштабированием, аналогично тому, как это было сделано для полиномов Чебышёва в §2.3б.

136 2. Численные методы анализа Помимо полиномов Лежандра существуют и другие семейства орто гональных полиномов, широко используемые в теории и практических вычислениях. В частности, введённые в §2.3 полиномы Чебышёва обра зуют семейство полиномов, ортогональных на интервале [1, 1] с весом (1 x2 )1/2.

Часто возникает необходимость воспользоваться ортогональными полиномами на бесконечных интервалах [0, +] или даже [, ].

Естественно, единичный вес (x) = 1 тут малопригоден, так как с ним интегралы по бесконечным интервалам окажутся, по большей части, расходящимися. Полиномы, ортогональные на интервалах [0, +] или [, ] с быстроубывающими весами ex и ex называются поли номами Лагерра и полиномами Эрмита соответственно.10 Они также находят многообразные применения в задачах приближения, и более подробные сведения на эту тему читатель может почерпнуть в [26, 62].

2.12 Численное интегрирование 2.12а Постановка и обсуждение задачи Задача вычисления определённого интеграла b (2.109) f (x) dx a является одной из важнейших математических задач, к которой сво дится большое количество различных прикладных вопросов. Это на хождение площадей криволинейных фигур, центров тяжести и момен тов инерции тел, работы переменной силы и т. п. механические, фи зические, экономические и другие задачи. В математическом анализе обосновывается формула Ньютона-Лейбница b (2.110) f (x) dx = F (b) F (a), a первообразная для функции f (x), т. е. F (x) = f (x). Она где F (x) даёт удобный способ вычисления интегралов, который в значительной степени удовлетворяет потребности практики. Тем не менее, возникают 10 Иногда их называют также полиномами Чебышёва-Лагерра и Чебышёва Эрмита (см., к примеру, [39, 62]), поскольку они были известны ещё П.Л.Чебышёву.

2.12. Численное интегрирование ситуации, когда для вычисления интеграла (2.109) требуются другие подходы.

Задачей численного интегрирования называют задачу нахождения определённого интеграла (2.109) на основе знания значений функции f (x), без привлечения её первообразных и формулы Ньютона-Лейбница (2.110). Подобная задача нередко возникает на практике, например, ес ли подинтегральная функция f (x) задана таблично, т. е. своими значе ниями в дискретном наборе точек, а не аналитической формулой. В некоторых случаях численное нахождение интеграла приходится вы полнять потому, что первообразная для интегрируемой функции не выражается через элементарные функции. Но даже если эта первооб разная может быть найдена в конечном виде, её вычисление не всегда осуществляется просто (длинное и неустойчивое к ошибкам округления выражение и т. п.). Все эти причины вызывают необходимость разви тия численных методов для нахождения определённых интегралов.

?

x Рис. 2.20. Вычисление определённого интеграла необходимо при нахождении площадей фигур с криволинейными границами Наибольшее распространение в вычислительной практике получили формулы вида n b (2.111) f (x) dx ck f (xk ), a k= где ck некоторые постоянные коэффициенты, xk точки из интерва ла интегрирования [a, b], k = 0, 1,..., n. Подобные формулы называют квадратурными формулами,11 коэффициенты ck это весовые коэфи циенты или просто вес квадратурной формулы, а точки xk её узлы.

а 11 Что буквально означает формулы вычисления площади.

138 2. Численные методы анализа В многомерном случае аналогичные приближённые равенства n f (x) dx ck f (xk ), k= где x, xk Rm, область в Rm, m 2, называют кубатурными формулами. Естественное условие принадлеж ности узлов xk области интегрирования вызвано тем, что за её предела ми подинтегральная функция может быть просто не определена, как, например, arcsin x вне интервала [1, 1].

Тот факт, что квадратурные и кубатурные формулы являются ли нейными выражениями от значений интегрируемой функции в узлах, объясняется линейным характером зависимости самого интеграла от подинтегральной функции. С другой стороны, квадратурные форму лы можно рассматриваеть как обобщения интегральных сумм Римана (через которые интеграл Римана и определяется), так как простейшие составные квадратурные формулы прямоугольников просто совпадают с этими интегральными суммами.

Как и ранее, совокупность узлов x0, x1,..., xn квадратурной (ку батурной) формулы называют сеткой. Разность n b f (x) dx R(f ) = ck f (xk ) a k= называется погрешностью квадратурной формулы или её остаточным членом. Это число, зависящее от подинтегральной функции f, в отли чие от остаточного члена интерполяции, который является ещё функ цией точки.

Если для некоторой функции f или же для целого класса функций F f имеет место точное равенство n b f (x) dx = ck f (xk ), a k= то будем говорить, что квадратурная формула точна (является точ ной) на f или для целого класса функций F. То, насколько широким является класс функций, на котором точна рассматриваемая формула, может служить косвенным признаком её точности вообще. Очень ча сто в качестве класса пробных функций F, для которых исследуется 2.12. Численное интегрирование совпадение результата квадратурной формулы и искомого интеграла, берут алгебраические полиномы. В этой связи полезно Определение 2.12.1 Алгебраической степенью точности квадратур ной формулы называют наибольшую степень алгебраических полино мов, для которых эта квадратурная формула является точной.

Сответственно, из двух квадратурных формул более предпочти тельной будем считать ту, которая имеет бльшую алгебраическую сте о пень точности. Неформальным обоснованием этого критерия служит тот факт, что с помощью полиномов более высокой степени можно по лучать более точные приближения функций, как локально (с помощью формулы Тейлора), так и глобально (к примеру, с помощью разложе ния по полиномам Чебышёва или Лежандра).

Рассмотрим теперь влияние погрешностей реальных вычислений на ответ, получаемый с помощью квадратурных формул. Предположим, что значения f (xk ) интегрируемой функции в узлах xk вычисляются неточно, с погрешностями k. Тогда при вычислениях по квадратурной формуле получим n n n ck f (xk ) + k = ck f (xk ) + ck k.

k=0 k=0 k= Если для всех k = 0, 1,..., n знаки погрешностей k совпадают со зна ками весов ck, то общая абсолютная погрешность результата, получен ного по квадратурной формуле, становится равной k |ck |k. Следова тельно, сумму модулей весов квадратурной формулы, т. е. величину n |ck |, k= нужно рассматривать как коэффициент усиления погрешности при вы числениях с этой формулой.

Условие устойчивости вычислений по квадратурной формуле при значительном количестве узлов n вызывает требование положительно сти весов ck, поскольку тогда при прочих равных условиях сумма мо дулей весов минимальна. В частности, в случае интегрирования функ ций, принимающих значения одного знака, мы избегаем тогда потери точности при вычитании близких значений.

140 2. Численные методы анализа 2.12б Формулы Ньютона-Котеса.

Простейшие квадратурные формулы Простейший приём построения квадратурных формул замена под интегральной функции f (x) на интервале интегрирования [a, b] на бо лее простую, легче интегрируемую функцию, которая интерполирует или приближает f (x) по заданным узлам x0, x1,..., xn. В случае, когда f (x) заменяется интерполянтом и все рассматриваемые узлы простые, говорят о квадратурных формулах интерполяционного типа, или, что равносильно, об интерполяционных квадратурных формулах.

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

Формулами Ньютона-Котеса называют интерполяционные квад ратурные формулы, полученные с помощью алгебраической интерпо ляции подинтегральной функции на равномерной сетке с простыми уз лами. В зависимости от того, включаются ли концы интервала инте грирования [a, b] в множество узлов квадратурной формулы, различа ют формулы Ньютона-Котеса замкнутого типа и открытого типа.

Построим формулы Ньютона-Котеса для n = 0, 1, 2, причём будем строить наиболее популярные формулы замкнутого типа за исключе нием случая n = 0, когда имеем всего один узел и замкнутая формула просто невозможна.

Если n = 0, то подинтегральная функция f (x) интерполируется по линомом нулевой степени, т. е. какой-то константой, равной значению f (x) в единственном узле x0. Если взять x0 = a, то получается квадра турная формула левых прямоугольников, а если x0 = b формула правых прямоугольников (см. Рис. 2.21).

Ещё один естественный вариант выбора единственного узла x0 = 2 (a + b), т. е. середина интервала интегрирования [a, b]. При этом приходим к квадратурной формуле b a+b f (x) dx (b a) · f, a называемой формулой средних прямоугольников: согласно ей интеграл берётся равным площади прямоугольника с основанием (b a) и вы 2.12. Численное интегрирование Рис. 2.21. Иллюстрация квадратурных формул левых и правых прямоугольников сотой f ((a + b)/2) (см. Рис. 2.22). Эту формулу часто называют так же просто формулой прямоугольников, так как она является наибо лее часто используемым вариантом рассмотренных простейших квад ратурных формул.

Оценим погрешность формулы средних прямоугольников методом локальных разложений, который ранее был использован при исследо вании численного дифференцирования. Разлагая f (x) в окрестности точки x0 = 2 (a + b) по формуле Тейлора с точностью до членов перво го порядка, получим f () a+b a+b a+b a+b + f · x · x f (x) = f +, 2 2 2 2 где зависящая от x точка интервала [a, b], которую корректно обо значить через (x). Далее b a+b f (x) dx (b a) · f R(f ) = a b a+b f (x) f = dx a b f ((x)) a+b a+b a+b f · x · x = + dx 2 2 2 a b f ((x)) a+b · x = dx, 2 a 142 2. Численные методы анализа Рис. 2.22. Иллюстрация квадратурных формул средних прямоугольников и трапеций поскольку ba b a+b x dx = t dt = 0, 2 ba a интеграл от первого члена разложения зануляется. Следовательно, с учётом принятого нами ранее обозначения Mp := max f (p) (x) x[a,b] можно выписать оценку b b f () 2 a+b M2 a+b |R(f )| · x dx x dx 2 2 2 a a 3b M2 (b a) M2 1 a+b · x = =.

23 2 a Отсюда, в частности, следует, что для полиномов степени не выше формула (средних) прямоугольников даёт точное значение интеграла, коль скоро вторая производная подинтегральной функции тогда зану ляется и M2 = 0.

Полученная оценка точности неулучшаема, так как достигается на функции g(x) = x 2 (a + b). При этом a+b M2 = max |g (x)| = 2, g = 0, x[a,b] 2.12. Численное интегрирование и потому b (b a)3 M2 (b a) a+b g(x) dx (b a) · g = =, 2 12 a т. е. имеем точное равенство на погрешность.

Рассмотрим теперь квадратурную формулу Ньютона-Котеса, соот ветствующую случаю n = 1, когда подинтегральная функция прибли жается интерполяционным полиномом первой степени. Построим его по узлам x0 = a и x1 = b:

xb xa P1 (x) = f (a) + f (b).

ab ba Интегрируя это равенство, получим b b b f (a) f (b) (x b) dx + (x a) dx P1 (x) dx = ab ba a a a b b f (a) (x b)2 f (b) (x a) = + ab ba 2 a a ba = f (a) + f (b).

Мы вывели квадратурную формулу трапеций b f (a) + f (b) (2.112) f (x) dx (b a) ·, a название которой также навеяно геометрическим образом. Фактиче ски, согласно этой формуле точное значение интеграла заменяется на значение площади трапеции (стоящей боком на оси абсцисс) с высотой (b a) и основаниями, равными f (a) и f (b) (см. Рис. 2.22).

Чтобы найти погрешность формулы трапеций, вспомним оценку (2.24) для погрешности интерполяционного полинома. Из неё следует, что f ((x)) f (x) P1 (x) = · (x a)(x b) для некоторой точки (x) [a, b]. Таким образом, для формулы трапе ций b b f ((x)) f (x) P1 (x) dx = · (x a)(x b) dx, R(f ) = a a 144 2. Численные методы анализа но вычисление полученного интеграла на практике нереально из-за неизвестного вида (x). Как обычно, имеет смысл вывести какие-то бо лее удобные оценки погрешности, которые, возможно, будут не столь точны.

Поскольку выражение (x a)(x b) почти всюду на интервале [a, b] сохраняет один и тот же знак, то b f ((x)) |R(f )| · |(x a)(x b)| dx a b M · (x a)(x b) dx, 2 a где M2 = maxx[a,b] |f (x)|. Далее b b x2 (a + b) x + ab dx (x a)(x b) dx = a a b b x3 x2 b (2.113) (a + b) = + abx 3 2 a a a 2 b3 a3 3(a + b) b2 a2 + 6ab(b a) = (b a) b3 + 3ab2 3a2 b + a3 = =.

6 Поэтому окончательно M2 (b a) |R(f )|.

Эта оценка погрешности квадратурной формулы трапеций неулучша ема, поскольку достигается на функции g(x) = (x a)2.

2.12в Квадратурная формула Симпсона Построим квадратурную формулу Ньютона-Котеса для n = 2, т. е.

для трёх равномерно расположенных узлов a+b x0 = a, x1 =, x2 = b 2.12. Численное интегрирование из интервала интегрирования [a, b].

Для упрощения рассуждений выполним параллельный перенос кри волинейной трапеции, площадь которой мы находим с помощью инте грирования, и сделаем началом координат оси абсцисс точку a. Тогда правым концом интервала интегрирования станет l = b a. Пусть P2 (x) = c0 + c1 x + c2 x полином второй степени, интерполирующий сдвинутую подинтег ральную функцию по узлам 0, l/2 и l. Если график P2 (x) проходит через точки плоскости с координатами l a+b 0, f (a),,f, l, f (b), 2 то c0 = f (a), l l a+b (2.114) c0 + c1 + c2 = f, 2 4 c0 + c1 l + c2 l2 = f (b).

С другой стороны, площадь, ограниченная графиком интерполяцион ного полинома, равна l l2 l c0 + c1 x + c2 x2 dx = c0 l + c1 + c 2 l 6c0 + 3c1 l + 2c2 l2.

= Фактически, для построения квадратурной формулы требуется решить относительно c0, c1 и c2 систему уравнений (2.114) и потом подставить результаты в полученное выше выражение. Но можно выразить 6c0 + 3c1 l + 2c2 l2 через значения подинтегральной функции f в узлах, не решая систему (2.114) явно.

Умножая второе уравнение системы (2.114) на 4 и складывая с пер вым и третьим уравнением, получим a+b 6c0 + 3c1 l + 2c2 l2 = f (a) + 4f + f (b).

146 2. Численные методы анализа a b 0 l Рис. 2.23. Иллюстрация вывода квадратурной формулы Симпсона Таким образом, l l c0 + c1 x + c2 x2 dx P2 (x) dx = 0 l a+b = f (a) + 4f + f (b), 6 что даёт приближённое равенство b ba a+b (2.115) f (x) dx · f (a) + 4f + f (b).

6 a Оно называется квадратурной формулой Симпсона или формулой па рабол (см. Рис. 2.23), коль скоро основано на приближении подинте гральной функции подходящей параболой. Читатель может самостоятельно убедиться, что та же самая форму ла получается в результате интегрирования по [a, b] интерполяционного 12 Приведённый нами элегантный вывод формулы Симпсона восходит к учебнику А.Н. Крылова [16].

2.12. Численное интегрирование полинома второй степени в форме Лагранжа a+b x (x b) (x a)(x b) a+b P2 (x) = f (a) + f a+b a+b a+b a (a b) a b 2 2 a+b (x a) x + f (b) a+b (b a) b 2 a+b a+b x (x b) f (a) 2(x a)(x b) f = (b a)2 2 a+b + (x a) x f (b), который строится для подинтегральной функции по узлам a, (a + b)/ и b.

Предложение 2.12.1 Квадратурная формула Симпсона имеет ал гебраическую степень точности 3, т. е. является точной для любого полинома степени не выше третьей.

Доказательство. Отметим прежде всего, что для полиномов степени не выше второй этот факт следует прямо из того, что формула строи лась нами как интерполяционная квадратурная формула, основанная на интерполяции подинтегральной функции полиномом второй степе ни. Поэтому достаточно показать, что формула Симпсона точна для монома x3, но не является точной для более высоких степеней.

Имеем b b 4 a x3 dx =.

a 148 2. Численные методы анализа С другой стороны, b a 3 a3 + 3a2 b + 3ab2 + b ba 3 a+b + b3 + b a +4 = a+ 6 2 6 b a 3a3 + 3a2 b + 3ab2 + 3b · = 6 b 4 a ba a + a2 b + ab2 + b3 = =, 4 что совпадает с результатом точного интегрирования.

Для монома x4 длинными, но несложными выкладками нетрудно проверить, что результат, даваемый формулой Симпсона для интегра ла по интервалу [a, b], т. е.

ba 4 a+b + b4, a + 6 отличается от точного значения интеграла b b 5 a x4 dx = a на величину (b a) /120. Она не зануляется при a = b, так что на полиномах четвёртой степени формула Симпсона уже не точна.

Итак, несмотря на то, что в формуле Симпсона интерполяционный полином для подинтегральной функции строился нами по трём точ кам и имел степень 2, фактическая точность формулы более высока, чем та, что обеспечивается полиномом второй степени. Для правильной оценки погрешности формулы Симпсона в этой ситуации нужно взять для подинтегральной функции интерполяционный полином более вы сокой степени третьей. При наличии всего 3 узлов мы находимся в условиях задачи интерполяции с кратными узлами.

Предполагая существование производной f в среднем узле x1 = (a + b)/2, можно считать, к примеру, что именно он является крат ным узлом. При этом мы будем формально решать задачу построения такого интерполяционного полинома 3-й степени H3 (x), что (2.116) H3 (a) = f (a), H3 (b) = f (b), a+b a+b a+b a+b = f (2.117) H3 =f, H3.

2 2 2 2.12. Численное интегрирование Конкретное значение производной в средней точке (a+b)/2 далее никак не будет использоваться. Здесь нам важно лишь то, что при любом значении этой производной решение задачи (2.116)–(2.117) существует, и потребуется оценка его погрешности.

Существование и единственность решения подобных задач мы уста новили в §2.4 и там же привели оценку его погрешности (2.44):



Pages:     | 1 | 2 || 4 | 5 |   ...   | 10 |
 





 
© 2013 www.libed.ru - «Бесплатная библиотека научно-практических конференций»

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