Численные методы. Соснин (2005) (1160462), страница 14
Текст из файла (страница 14)
Из принятых условий следует,что в точке λ = 0 функция g(λ) достигает минимума. Она, очевидно, дифференцируема, поэтомуg 0 (0) = 0. Продифференцировав (4.20) по λ, получим, чтоAc − f , y = 0.2λ hAy, yi + 2 y, Ac − f λ=0 = 0 ⇐⇒Из произвольности выбора y следует, что Ac − f = 0. Это означает, что c — решение системы(4.19). Необходимость доказана.Осталось заметить, что из доказанной эквивалентности следует существование и единственностьточки минимума функции F (c). Это вытекает из того, что матрица A положительно определена,поэтому система (4.19) имеет единственное решение. Теорема полностью доказана.Алгоритм построения наилучшего приближенияОпираясь на доказанную теорему, можно построить алгоритм нахождения наилучшего приближения для функции f ∈ L2 [a; b].
Он будет выглядеть так:1. Выбираем (n + 1) линейно независимый элемент ϕk , k = 0, n из L2 [a; b].2. Строим матрицу A = (akl ) скалярных произведений:Zbakl =ϕk (x)ϕl (x) dx.a3. Готовим вектор скалярных произведений f = (f0 , f1 , . . . , fn ), где fi находятся так:Zbfi =f (x)ϕi (x) dx,i = 0, n.a4. Ищем вектор коэффициентов c = (ec0 , ec1 , . . .
, ecn ), решая систему уравнений Ac = f .5. Строим элемент ϕe:ϕe=ec0 ϕ0 + ec1 ϕ1 + . . . + ecn ϕn .Он будет являться наилучшим приближением согласно доказанной теореме.Теперь посмотрим, насколько точно ϕe приближает f, то есть оценим ||f − ϕ||.e Для этого нампонадобится лемма.Лемма. Пусть ϕe — элемент наилучшего приближения для f. Тогдаhf − ϕ,e ϕie =0— то есть ϕe ортогонален (f − ϕ).eДоказательство. Подставим представление ϕe через ϕk в искомое скалярное произведение:*+nnnn XnXXXXhf − ϕ,e ϕie = f−eck ϕk ,ecl ϕl =ecl hf, ϕl i −eck ecl hϕk , ϕl i =k=0=nXecl fl −l=0n XnXl=0l=0k=0 l=0eck ecl akl = {согласно старым обозначениям } = c, f − hAc, ci = f − Ac, c .k=0 l=0Согласно построению ϕ,e c находился из условия Ac = f , поэтому получаем, чтоhf − ϕ,e ϕie = 0.Лемма доказана.4.2.
НАИЛУЧШЕЕ ПРИБЛИЖЕНИЕ ТАБЛИЧНОЙ ФУНКЦИИ71Теперь можно оценить отклонение:||f − ϕ||e 2 = hf − ϕ,e f − ϕie = hf − ϕ,e f i − hf − ϕ,e ϕie = hf − ϕ,e f i = hf, f i − hϕ,e fi .Согласно лемме, hf, ϕie = hϕ,e ϕie , поэтому||f − ϕ||e 2 = hf, f i − hϕ,e ϕie = ||f ||2 − ||ϕ||e 2.Замечание. Если {ϕk } — ортонормированная система, то есть hϕek , ϕel i = δkl , тогда A = E. Отсюдаследует, что ck = hf, ϕk i = fk , и для наилучшего приближения получается простая формула:ϕe=nXfk ϕk .i=0В этом случае коэффициенты ck называются коэффициентами Фурье, а построенный элементϕe — многочленом Фурье.Пример 4.3.Построим для функции из предыдущего примера наилучшее приближение.
Итак, f (x) задана таблично в точках x0 , x1 = x0 + h, x2 = x0 + 2h. Обозначим F0 = f (x0 ), F1 = f (x1 ), F2 = f (x2 ).Приближать будем снова многочленами:ϕ0 (x) = 1,ϕ1 (x) = x − x1 .Напомним, что, применяя метод наименьших квадратов, мы нашли такое приближение:Φ(x) =F2 − F0F0 + F1 + F2+(x − x1 )32hДля поиска наилучшего приближения для f последовательно пройдем по построенному нами алгоритму:1.
ϕ0 (x) = 1,ϕ1 (x) = x − x1 .xZ1 +h2. Подсчитаем akl =ϕk (x)ϕl (x) dx :x1 −ha00=2h;xZ1 +ha10= a01 =x1 −hxZ1 +ha11(x − x1 )2 dx ==x1 −h2h=⇒ A = 0x +h(x − x1 )2 1= 0;(x − x1 ) dx =2x1 −hx +h(x − x1 )3 1h32h3h3+=.=3333x1 −h02h3 .33. На данном этапе возникают сложности, так как, не зная f (x), мы не можем точно вычислитьf0 и f1 . Будем вычислять их приближенно: f0 через формулу среднего значения, а f1 — по72Глава 4. ИНТЕРПОЛЯЦИЯ И ПРИБЛИЖЕНИЕ ФУНКЦИЙформуле Симпсона (как известно, она дает маленькую погрешность):xZ1 +hf0f (x) dx = f (x) · 2h ≈ 2h=F0 + F1 + F2;3x1 −hxZ1 +hf1Zbb−af (x)(x − x1 ) dx ≈ { G(x) dx ≈(G0 + 4G1 + G2 )} ≈6=ax1 −h≈h22h[F0 (−h) + 4F1 · 0 + F2 · h] =(F2 − F0 )634.
Теперь решаем систему Ac = f . Матрица A — диагональная, поэтому ее решение запишетсяпросто:F0 + F1 + F2F0 + F1 + F2c0 = 2h;c0 =; 2he e33⇐⇒ 2h3F2 − F0h2 ec1 =.ec1 =(F2 − F0 ).2h335. Вычислив коэффициенты ec0 и ec1 , можем записать построенное приближение:F0 + F1 + F2F2 − F0ϕ(x)e=+(x − x1 ).32hОно совпало с построенной ранее функцией Φ(x).
Совпадение это не случайно: мы очень неточновычислили f0 , хотя могли бы этого избежать. Вычислим f0 , применяя формулу Симпсона:xZ1 +hf (x) dx ≈h[F0 + 4F1 + F2 ].3x1 −hЗаново решив систему Ac = f , получим такие выражения для ec0 и ec1 :F+4F+F012c0 =; e6F2 − F0 ec1 =.2hЭти коэффициенты дадут более точное приближение:F0 + 4F1 + F2F2 − F0ϕ(x)e=+(x − x1 ).62hЕстественно задаться вопросом: «А насколько оно точнее?». Легко показать, что ϕe отличается отϕe на константу:F0 − 2F1 + F2h2 F0 − 2F1 + F2h2ϕ(x)e− ϕ(x)e==·=fxx,1 ,66h26где fxx,1 — вторая разностная производная f (x) в точке x1 .
Ранее было показано, что такое жезначение принимает квадрат нормы погрешности между ϕe и f — то есть ||ϕ(x)e− f (x)||2 :||ϕ(x)e− f (x)||2 =1(F0 − 2F1 + F2 )2 .6Неформально можно сказать, чтоf −ϕe≈ϕe − ϕ.eСтоль большие отклонения возникают из-за неточности при вычислении f . Отсюда следует вывод, что алгоритм построения наилучшего приближения слишком зависит от методов приближенныхвычислений, и лучше использовать единые формулы, например, для подсчета интегралов — ту жеформулу Симпсона. При этом получаются неплохие результаты.Глава 5Численные методы решения краевыхзадачНашей первой задачей будет поиск численных методов решения задачи Коши:du= f (t, u(t)), 0 < t < T ;dt u|= u .t=00Для начала необходимо построить соответствующую дискретную модель. Для этого разобьем весьотрезок [0; T ] на точки ωτ = {tn = nτ, n = 0, 1, .
. . Tτ }, где τ — диаметр дискретной сетки. Обозначимзначения искомой функции — un = u(tn ), приближенное решение — yn , и погрешность на n-й итерациикак zn = yn − un в узлах сетки (zn , yn , un — сеточные функции).От будущего алгоритма потребуем как можно более точного воспроизведения функции u — дляэтого нужно, чтобы погрешность zn была мала.Обсудим понятие сходимости приближенного решения к точному. Фиксируем точку tn и построимпоследовательность сеток ωτ такую, что точка tn является узлом для сеток с номерами m > k, тоесть при сгущении сетки только добавляются новые узлы. На каждой из этих сеток строится сеточнаяфункция yn .τ →0Определение.
Сеточная функция yn сходится к решению un в узле tn , если |zn | −→ 0.Определение. Сходимость на отрезке означает сходимость в каждой точке этого отрезка.Определение. Пусть погрешность по порядку роста ведет себя как |zn | = O(τ p ), тогда приближенное решение имеет p-й порядок точности.Если мы будем приближать производную в узлах сетки ее разностным аналогом:yn+1 − ynu0 (tn ) =,τто исходное уравнение примет такой вид:yn+1 − yn= f (tn , yn ), n = 0, 1, . . .(5.1)τПолучившиеся уравнения относительно yi называются разностной схемой.Из этого аналога нашего непрерывного уравнения можно определить значения yn во всех точкахсетки, если стартовать с фигурирующего в условии задачи y0 .yn+1 = yn + τ f (tn , yn ); n = 0, 1, .
. .y0 = u 0 .7374Глава 5. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ КРАЕВЫХ ЗАДАЧПолучившаяся схема называется явной, так как (n+1)-е приближение выражается через n-е явно,непосредственно (не требуется решать уравнение).Рассмотрим другой вариант, в котором мы заменяем обыкновенную производную на разностнуюпроизводную назад. Тогда функция f (t, y) будет браться в точке (tn+1 , yn+1 ) :yn+1 − yn= f (tn+1 , yn+1 ).τПри этом относительно yn+1 возникнет нелинейное уравнение:yn+1 = τ f (tn+1 , yn+1 ) + yn , y = 0, 1, . . . ;y0 = u 0 .То есть на каждом шаге требуется решать нелинейное, вообще говоря, уравнение, и формула дляy соответственно является неявной. Отсюда и название — неявные разностные схемы.un+1 − un+ f (tn , un ) называется невязкой (погрешностьюОпределение.
Значение ψn = −τаппроксимации) расчетной схемы (5.1) в узле tn .Примечание. Вообще говоря, далее под невязкой мы будем понимать разность левой и правой частей расчетной схемы при подстановке в нее точного решения. Вышеприведенная формула выражает невязку для методов Рунге-Кутта с параметром m = 1.τ →0Определение.
Разностная схема аппроксимирует исходное ОДУ в узле tn , если ψn −→ 0 (этотпредельный переход связан с последовательностью сеток).Определение. Разностная схема имеет p-й порядок аппроксимации в точке tn , если выполненоравенство: |ψn | = O(τ p ).Данную терминологию будем использовать для сравнения методов.5.1Сходимость методов Рунге-КуттаМетодами Рунге-Кутта1 называется семейство методов, общий вид разностных схем которых задаетсятак:mXyn+1 − yn=σi Ki (y),(5.2)τi=1где величины Ki вычисляются по следующим формулам:K1 (y)=f (tn , yn );K2 (y)=f (tn + a2 τ, yn + b21 τ K1 (y));K3 (y)= f (tn + a3 τ, yn + b31 τ K1 (y) + b32 τ K2 (y));...m−1XKm (y) = f (tn + am τ, yn +bmi τ Ki (y)).i=1Параметры ai , bij , σi выделяют конкретный метод Рунге-Кутта.В этих методах нужно вычислять значения функции f в промежуточных точках сетки.
Их количество определяется параметром m, а соответствующие методы называются m-этапными методами.Схемы с m > 5 используются крайне редко (чаще всего используют 4-этапные методы).1 Основнаяидея — К. Рунге (1885).5.1. СХОДИМОСТЬ МЕТОДОВ РУНГЕ-КУТТА75Какие ограничения накладываются на параметры методов Рунге-Кутта для того, чтобы обеспечитьсходимость? Попробуем ответить на этот вопрос.Потребуем, чтобы разностная схема (5.2) аппроксимировала исходное ОДУ в соответствии с ввеτ →0денным определением ( |ψn | −→ 0).
Невязка в данном случае будет иметь вид:mψn = −un+1 − un X+σi Ki (u).τi=1Разложим un+1 в ряд Тейлора с центром в точке tn :un+1 = un + u0n τ + O(τ 2 ).Рассмотрим выражение для Ki (аналогично разложив его в ряд Тейлора):Ki (u) = f (tn , un ) + O(τ ).Подставим эти формулы в выражение для невязки:ψn =−u0n+ f (tn , un )mXσi + O(τ ) = f (tn , un ) −1 +i=1mX!σi+ O(τ ).i=1Очевидно, для того, чтобы разностная схема аппроксимировала исходное ОДУ (с порядком апmXσi = 1; тогда невязка будет равнапроксимации p = 1, то есть ψn = O(τ )), достаточно, чтобыi=1O(τ ).Таким образом, можно ввести первое ограничение на параметры метода (5.2):mXσi = 1.i=1Теорема 5.1 (О сходимости методов Рунге-Кутта). Пусть метод Рунге-Кутта аппроксимируетисходное уравнение, тогда приближенное решение yn сходится к точному un , и порядок точностиприближенного решения совпадает с порядком аппроксимации разностной схемы исходного ОДУ.Доказательство.