Г.М. Кобельков - Курс лекций по численным методам (1160467), страница 6
Текст из файла (страница 6)
А этои означает, что функция является линейной.Для более высоких степеней проще использовать соответствующие результаты вариационного исчисления.Для экстремалей функционалаZbΦ(x) := F (t, x, ẋ, . . . , x(k) )dt → min(82)aимеет место уравнение Эйлера – Пуассона:ddddFẋ −Fẍ − · · · −Fx(k−1) − Fx(k) . . .= 0.Fx −dtdtdtdt(83)В частном случае, при k = 2, оно превращается в уравнениеd2dFẍ − Fẋ + Fx = 0.2dtdt(84)В нашем случае F = (ẍ)2 , поэтому экстремаль функционала I2 удовлетворяет уравнениюd2 ′′s (x) = 0,dx2(85)то есть получаем условие экстремума s(4)(x) = 0. Это значит, что экстремаль — многочлен степени не более 3.Введём обозначение hn+1 := xn+1 − xn . Следовательно, s′′ (x) имеет видs′′ (x) = Mnxn+1 − xx − xn+ Mn+1.hn+1hn+119(86)У нас есть условия s(xn ) = f (xn ) и s(xn+1 ) = f (xn+1 ).
Можно показать, что в этом случае s имеет видs(x) =Mn (xn+1 − x)3Mn+1 (x − xn )3xn+1 − xx − xn+·+ An+ Bn.6hn+16hn+1hn+1hn+1(87)Подставляя узлы x = xn и x = xn+1 , получаем уравнения на коэффициенты:An = f (xn ) −Mn 2h,6 n+1Bn = f (xn+1 ) −Mn+1 2hn+1 .6(88)Осталось найти Mn , и его мы найдём из соображений непрерывности первой производной. Напишем формулудля производной сплайна для отрезков [xn−1 , xn ] и [xn , xn+1 ]:s′ (x) = −Mn−1 (xn − x)2Mn (x − xn−1 )2An−1Bn−1+−+,2hn2hnhnhn(89)Mn (xn+1 − x)2Mn+1 (x − xn )2AnBn+−+,2hn+12hn+1hn+1hn+1(90)s′ (x) = −Затем подставим в оба равенства узел xn :1 Mn 2Mn−1 2hn − f (xn−1 ) +hn + f (xn ) .s′ (xn ) =hn361Mn 2Mn+1 2s′ (xn ) =−hn+1 − f (xn ) −hn+1 + f (xn+1 ) .hn+136(91)(92)Итого получим такое уравнение:Mn+1hn+1 + hnMn−1f (xn+1 ) − f (xn ) f (xn ) − f (xn−1 )hn+1 + Mn+hn =−= f (xn+1 ; xn ) − f (xn−1 ; xn ).
(93)636hn+1hnПолучаем систему с трёхдиагональной матрицей, у которой имеется диагональное преобладание. В самом деле,hn+1 + hn >hn+1hn+.22(94)Впрочем, пока здесь есть одна проблема. Это равенство справедливо только в промежуточных узлах, и поэтому нам не удастся получить замкнутую систему линейных уравнений. Для замыкания системы нужно ещёдва уравнения, которые приходится брать «с потолка».
Например, можно сделать так: провести интерполяционный многочлен Лагранжа L через узлы x0 , x1 , x2 , x3 и положить M0 := L′′ (x0 ), потому что коэффициенты Miв каком-то смысле аппроксимируют значения производной исходной функции. Аналогично поступаем и с точкой xn .Определение. Сплайны такого типа называются интерполяционными.Для сплайнов второго порядка имеют место оценки: (k)f − s(k) 6 Ch4−k , k = 0, .
. . , 3, C = f (4) , h = max hi ,(95)то есть приближается не только сама функция, но и её производные.Но сразу видно, чем эти сплайны плохи. При добавлении всего одной точки в систему узлов придётся всекоэффициенты пересчитывать заново, что очень неудобно.
Значит, нужно придумывать что-то более локальное.О нём пойдёт речь в следующем параграфе.2.7.2. В-сплайнЭтот сплайн изобретён как попытка бороться с недостатками интерполяционных сплайнов. Правда, приходится отказаться от того, чтобы сплайн проходил через точки (xi , f (xi )).Дальнейшая конструкция не является очень естественной, поэтому мы приведём только готовый рецепт, непоясняя, как до этого додумались.Рассмотрим отрезок [−2, 2] и разделим его на маленькие отрезки длины 1.
Рассмотрим функцию B(x),удовлетворяющую таким свойствам:• supp B = [−2, 2];• на каждом маленьком отрезке B является кубическим многочленом;• B ∈ C2 (R);20• B(x) чётна;• выполнено нормирующее условие B(0) + 2B(1) = 1.Если немного подумать, можно придумать, например, такую функцию3122 3 − x + 2 |x| , |x| 6 1;B(x) := 16 (2 − |x|)3 ,|x| ∈ [1, 2];0,|x| > 2.(96)А теперь построим кубический сплайн с использованием этой функции.Рассмотрим разбиение отрезка [a, b] на N частей с постоянным шагом h. Рассмотрим функцию(i)B3 (x):=N+1Xαn(i) Bn=−1x − nhh(97).(i)Это и будет наш сплайн. Пока мы не определили, что такое αn . Существуют два вида сплайнов.Первая разновидность (i = 1). Нам заданы значения функции f0 , .
. . , fN в узлах. Тогда с помощью линейной(1)интерполяции строятся значения f−1 и fN +1 . Затем полагаем αn := fn и рисуем наш сплайн.Вторая разновидность (i = 2). Для этого мы сначала с помощью кубической экстраполяции определяем значения функции f−2 , f−1 и fN +1 , fN +2 , проводя кубические сплайны через наборы точек x0 , . . . , x3 и xN −3 , . .
. , xN(2)соответственно. Затем определяем коэффициенты αn по формуламα(2)n :=8fn − fn+1 − fn−1.6(98)Почему нужно брать именно такие, а не другие, тайна сия велика есть.Но как бы там ни было, для этих сплайнов имеют место оценки: (k)f − (B (1) )(k) 6 Ch2−k , k = 0, 1;3 (k)f − (B (2) )(k) 6 Ch4−k ,3(99)(100)k0, . . . , 3.Основное преимущество этих сплайнов таково: они локальны.
То есть если значение функции меняется в однойточке, то не нужно пересчитывать все коэффициенты.3. Численные методы и дифференциальное исчисление3.1. Численное дифференцированиеБудем считать, что заданы узлы x1 , . . . , xn и значения функции f (x1 ), . . . , f (xn ). Самое простое — это свестизадачу к предыдущей, то есть приблизить функцию интерполяционным многочленом Лагранжа и считать, чтоL′n (x) ≈ f ′ (x).
Однако пример Адамара fε (x) = f (x) + ε sin εx2 показывает, что при стремлении к нулю в нормеC погрешности, норма производной погрешности стремится к бесконечности. Отсюда следует, что оператордифференцирования плохой и, вообще говоря, неограниченный.Для борьбы с плохими, сильно осциллирующими функциями, можно применять фильтрацию частот. Именно,для фильтрации спектра, сосредоточенного на отрезке [a, b], можно к функции f применять преобразованиеFb−1 χ[a,b] Fb .Выведем формулы для приближения производных в виде линейных комбинаций значений функции в некоторых точках, то естьXf ′ (x) ≈ai f (xi ).(1)Рассмотрим формулу f ′ (x) =f (x+h)−f (x)h+ r(x, h). Имеем r(x, h) =f ′′ (ξ)2· h, поэтому |r| 6hmax |f ′′ (x)|.2 [x,x+h]Таким образом, эта формула приближает производную с точностью до O(h).(x−h)Можно также применять формулу f ′ (x) ≈ f (x+h)−f.
Её точность — O(h2 ).2hПокажем, как строить формулы для аппроксимации производных по более чем двум узлам. Построим, например, формулу для первой производной по 3 узлам:f ′ (x) ≈1a1 f (x) + a2 f (x + h) + a3 f (x + 2h) .h21(2)Будем подбирать коэффициенты такими, чтобы формула давала точное значение производной для многочленовстепеней не выше 2. Получаем систему уравнений:a1 + a2 + a3 = 0,(3)a2 + 2a3 = 1,a2 + 4a4 = 0.Её решение — a1 = − 32 , a2 = 2, a3 = − 12 .
Значит, формула имеет видf ′ (x) ≈1−3f (x) + 4f (x + h) − f (x + 2h) .2h(4)Аналогично получим формулу для второй производной:f ′′ (x) ≈f (x + h) − 2f (x) + f (x − h) h2 (4)− f (x) + O(h4 ).h212(5)Оценим погрешность. Пусть |f ′′ (x)| 6 M , а ε — равномерная погрешность измерения функции, то есть |f −fe| 6 ε.2εТогда |f ′ − fe′ | 6 Mh2 + h . Чтобы погрешность измерений не сильно влияла на оценку производной, берёмpεh = 2 M . Отсюда√√√2M ε 2ε M∆6 √+ √ = 2 M ε.(6)2 ε2 M3.2. Сжатие информацииЗдесь была рассказана стандартная сказка про разложение по ОНБ в гильбертовом пространстве.
Основнаямысль: если хранить коэффициенты Фурье, то можно сильно сэкономить. Во-первых, часть коэффициентов(высокие частоты) можно выкинуть без особой потери смысла, ибо обычно они всё равно стремятся к нулю.Более того, набор коэффициентов можно более эффективно сжимать обычными методами — для не слишкомпоганых функций он выглядит куда приличнее, чем сама функция.3.2.1. Двумерный случайЕщё было что-то про сингулярное разложение матриц.
Какое оно имеет отношение к реальным алгоритмамсжатия, автор конспекта понимает с большим трудом. Точнее говоря, мы не будем сжимать никакие данные,а мы просто будем приближать функцию некоторыми базисными функциями, а хранить лишь коэффициентыразложения.Определение. Рассмотрим самосопряжённый оператор A. Разложение Σ = U t AV , где матрица Σ диагональна, а U и V ортогональны, называется сингулярным.Поставим такую задачу: найти наилучшее приближение для функции f ∈ L2 ([0, 1]2 ) с помощью функциивида ue(x)ev (y).
В качестве нормы в данном случае будем использовать L2 -норму. Положимu :=Как уже было сказано, мы хотим решить задачуue,keukv :=ve.kevk2(7)kf − µ · u(x)v(y)k → min .(8)ZZZ22u(x)f (x, y)v(y)dy dx + µu (x)dx · v 2 (y)dy.(9)Распишем эту норму:22kf − µuvk = kf k − 2µZВыражение в скобках под интегралом — это интегральный оператор с ядром f , применённый к функции v.Обозначим этот оператор через A.
Далее, заметим, что последнее слагаемое равно µ2 , потому что каждый изинтегралов равен 1 в силу выбранной нормировки функций u и v. С учётом всего этого можно переписать нормув таком виде:222kf − µ · uvk = kf k + µ − (Av, u) − (Av, u)2 .(10)Дальнейшая идея такова: будем максимизировать число (Av, u), чтобы вычитать побольше, а потом выберемтакое µ, чтобы второе слагаемое умерло, а именно, потом положим µ := (Av, u).22Реализуем этот план. Сначала оценим скалярное произведение:(11)|(Av, u)| 6 kAvk · kuk = kAvk .2Нам нужна функция v, для которой достигается максимум kAvk.
или, что то же самое, максимум kAvk . Имеем2kAvk = (Av, Av) = (A∗ Av, v).(12)Стало быть, если v отвечает максимальному собственному значению оператора A∗ A, то это число будет максимальным. Такой её и выберем.Остаётся вот какой вопрос: как найти это максимальное собственное значение. Для этого можно сделать сингулярное разложениематрицы, и выбрать максимальный по модулю диагональный элемент матрицы Σ.Алгоритм осуществления сингулярного разложения на лекциях рассказан не был.
Его можно сдуть из книжки Уилкинсона идр. «Справочник алгоритмов по линейной алгебре». — М., 1976.Фактически этот алгоритм эквивалентен полному QR-алгоритму нахождения собственных значений. Это ужасно долго. Еслинемного подумать и вспомнить доказательство одной теоремы из функана, то можно приближённо найти максимальное собственноезначение иным (итерационным) путём. А именно, берём любой единичный вектор, и начинаем к нему применять оператор A. Междуитерациями будем нормировать вектор, чтобы он не раздулся до бесконечности. Ежу ясно, что он будет стремиться к вектору смаксимальным собственным значением (попробуйте доказать или привести контрпример, возможно, получится и то и другое).
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
