Бабенко - Основы численного анализа (947491), страница 85
Текст из файла (страница 85)
Заметим, что не следует опасаться таких дискретизаций, которые приводят не к разреженным матрицам, а к полностью заполненным. Алгоритмы без насыщения, которые сами настраиваются на нужную 11окажем, что если т(0) б С (Я ), то и К(д, 0) б С,Я х Я ). В самом деле, К(д, 0) бесконечно дифференцируема по обеим переменным вне диагонали д = О. Далее, замечая, что т(д) = т(0) - (д — В)т (В) (~ у — ти(0) +..., ползучим 415 з 5, Крбатррвьге Формулы гладкость, позволяют получать решение более экономно, нежели алгоритмы невысокого порядка точности, но приводягцие к системам с разреженными матрицами, Примеры этому мы приведем ниже.
Качество дискретизации оценивают по свойствам матриц той линейной системы, к которой сводят задачу (речь идет о линейных задачах), Пусть Х= (о,,к(Оп О,)) и пусть 1 — единичная (2п-, 1) х (2п+ 1)-матрица. Важнейшей характеристикой матрицы Г з- М" . матрицы системы (24) — является ее мера обусловленности, С этой характеристикой мы подробнее ознакомимся в гл.
8, а сейчас отметим, что мера обусловленности произвольной квадратной матрицы А равна сопс(:1 = ~А~ А '~~, где ~ ~ — какая-либо норма матрицы. Найдем чебьппевскую норму матрицы .уе': Ос',оо = шах ~ ~~К(00 Оь)~ = шах — / )К(Оп О))с(0 < оо. 2 ! Г ь-о о Так как погрешность квадратурной формулы можно сделать сколь угодяо малой, то,1+Я'~ < Со Предгюженне, 3. Ьсли и > по, то ~(У+.Уб) '~ < Сг, где коистаигоа Сг, яе зависит от п. Доказатяльство. Прежде всего отметим, что ядро Л(0, 0) й С~;Б х хЯ'). В самом деле, нз уравнения (21) следует принадлежность ядра ЦО, О) как функции 0 пространству Сы(Я', Но в уравнении (2Ц можно поменять местами функции К и Ь, что равносильно тому, что мы будем рассматривать уравнение (20) как исходное, а уравнение (19) как резольвенту.
Тогда получим, что ядро 5(0, О) как функция О принадлежит пространству С Я'). Допустим, что матрица Г г.,Ж' особая; тогда найдется вектор к = (Хо,, Хг ) (,Х =- 1) такой, что (1г-,Х')ж = О. Построим по формуле (22) полинам 1(0) = 1(0, зс), и пусть Г(0) правая часть уравнения (19), в которозл вместо й мы подставили полинам 1(0).
Тогда Г(О) сумма полинома и гладкой функции. Беря ограничение полученного уравнения на сетку и используя то, что 1(0) — полипом, получим соотношение (23), в котором в правой части будет стоять Г(Вг). Заменяя интеграл суммой, получим уравнение вида зс + Хзс — —. 1 + р, где. у" =- Ц(Во),, Г(Вг,,)), р = О(п "), Отсюда Г =- О(п '). Обратимся к формуле (20) и возьмем ее ограничение на множество уъ~ов. Таь как Г сумма полинолга и гладкой функции, то -~ 5(0,0) Г(О)00 = о г г„ 2 —" — / Б(0 О) ~ Х(Оо)П (Π— Оь)ОВА. 0(п ') = 0(п ), о ь.=..о и, значит,к = 0(п "),что неверно. Следовательно,матрица Г + д~' пеособая. Глава О, Численное интегрирование Пусть вектор х = Яе,..., .Хзп) Щж = 1) таков, что ((Е+ Х) Х~ = ~(Е+Х) '/ . Построим паолинам с(0) = 1(0; х) и рассмотрим уравнение (19) с правой частью 1(0). Решение этого уравнения обозначим через ДО), и пусть й = (Га, ..., Гм,), Если в формуле (20) вместо Е подставить полянам ~, .то получим, что решение Г равно сумме полинома и гладкой функпииь Поэтому, делая дискретизацию уравнения (19), в котором правая часть заменена полиномом 1, получим (Е+Х)С = Х-га, где г = 0(тг ").
Отсюда й = (Е+Х) ~Х, (Е+Х) ~с, и поэтому > )(Е Ч Х ) '~ (1 — С . "). (25) С другой стороны, беря ограничение уравнения (20) на сетку и производя дискретизацию интеграла, получим ь = (1+.'2')Х вЂ”, О(п "), где .У вЂ”.— ( „',, ЦОн Оь)) . Поскольку с ь=е , ~)Е,(Ои О,)! = — ~~Ы(Он О)! 40, ь=о а то ~.х.'~ ( Са и, следовательно, ф ( 1 Р Са -- 0(п "'), откуда в силу оценки (25) следует утверждение предложения. Г) Итак, мы доказан|, что мера обустовленностн матрицы Š—; Х ограничена константой, не зависящей от п. Это очень важно, потому что при численном решении системы уравнений с данной матрнцей рост погрешности округления тем меньше, чека меньше мера обусловленности матрицы системы.
Подчеркнем, что в разностных методах решения краевых задач либо в методе конечных элементов потучаются системы с матрицами, мера обусловленности которых является величиной = 6 , где м ) 0 — некоторая постоянная. Этот факт сам по себе весьма неблагоприятный. Систему (24) можно решать либо методом Гаусса, если 2п -~- 1 невелико (2п, + 1 ( 121 + 141 для БЭСМ-6), либо методом простой итерации. Репшв эту систему, найдем плотность потенциала в узлах. Казалось бы, задача решена, поскольку остается только вычислить функпи|о (16), используя для этого формулу 1 /' Х(0)т'(0)гЕО к1 Е т(0) — в а я заменяя й(0) ицтерполяционным многочленом (22).
Но на самом деле основные трудности впереди. Вычислить этот интеграл без существенной потери точности, когда я близко к Ю, удается лишь после тщательной его регуляризации. '1то же касается вычисления производных эд(я), то положение еще хуже (по этим вопросам см. [49)). Мы же рассмотрим вопрос об определении функции а1 '(С), дающей отображение круга К на область Р. Вначале из уравнения (17), предварительно отделив мнимую часть и используя принятую параметризацию, найдем по заданному углу О точку - б ~7.
'!'огда получим 417 ' э5, Крбатррьые формулы уравнение Π—.— аг8(г(д) — яо) — — ~ х(5) Ве 05+ С. (2б) 1 0 т'(б) о Последний интеграл понимается как интеграл в смысле Коши. Рассмотрим вспомогательное ядро Л(д л) — Ве тЯ - т~(0) 218[(5 — д) ~2] Элементарно проверяется, что это ядро не имеет особенностей при 5 = д. Вне диагонали б = д оно бесконечно дифференцируемо по обеим переменным. Разлагая функцию т(5) в ряд Тейлора, нетрудно убедиться, что В(д, Я) б б С [Я' х Я'[. Поэтому соотношение (26) можно переписать в виде 0=:( (О)-:.) — -у1 у(о .г / 218[(б — д)/2~ 3 — — / ~~ — х(6 5+бб 1 / ( т'(б) х / ( тЯ вЂ” т(0) 2йй[(б — д)72) л подставляя вместо 1с(5) мпогочлсп (22), а затем, так жс как и выше, переходя От интеграла к сумме, получим 2 0= 8( (О) — ) — ', 1,Х,0 (О- 6,)-- а=о тЯ) 1 — +1 ~да[Во ( ) (О) — 2 [( ) 1 ~ ~ +С-~-0(п ").
(27) Вычислим предварительно константу С по формуле 1ш ~(ло) 4 С = О или по формуле о Затем в уравнении (27) отбросим величину 0(п ") - остаточный член аппроксимации, возпикгпей от замены интеграла на сумму. Полученное уравнение можно решать относительно д, задавая значения угла 0 = Оь = 2х1/(2п — , '1) (0 = О, 1,..., 2п). При решении уравнения применим метод Ньютона; проблем с определением начального приближения в данном случае не возникает. Найдя величины дь = д(Оь), находим соответствующие точки контура Уу: хь = = т(дь) — — (год)(0ь) Величины (зь) образуют таблицу функции я = (с) -- это ее значения в узлах бь = ехр(10ь1. Дая восстановления функции з =- з(Д воспользуемся формулой (5.3,43) (нужно только в этой формуле изменить обозначения). Для вычисления производной 'Я нужно продифференцировать формулу (22) (см.
также [187, 188, 191, 192, 195]). Тем саъ~ым конформное отображение круга па область построено. Этот метод изложен в [49), где представлены Глава б, Численное ннп!егрпроеапие многочиш!енные результаты расчетов. Сделаем два замечания практического порядка. Мы предположили, что контур области с? закан параметрическими уравнениями. Однако параметризация делается не единственным образом, поскольку после применения произвачьного диффеомо?эфизма отрезка !О, 2я) на себя мы снова получим некоторую параметризацию кривой Р. Удачная параметризация позволяет очень сильно повысить точность расчетов. Параметризацию высокого качества можно псшучить итерационным способом, решив предварительно задачу на небольшом числе узлов, а затем после анализа то!пь насколько хорошо решено уравнение (19), можно сделать вывод о том, как исправить параметризацию И второе замечание †.
как проверить точность расчетов, Один из тестов — это Таблица 1 опенка величины невязки (5.3.40). В данном случае удобно брать невязку не в видо (5,3.40), а следить за величинами 2 е! =- ~ ь ехр(г(йь), ь — — О 1 = 1, 2, ..., и, а точнее> определять величину е = твх !е!,'. ! Приведем данные одного расчета (табл. 1). Бралась область В, ограниченная кривой и = я+ аз~, где а константа,,а' < 1!!4, з = ехр(Ю), 0 < д <'?т. Функция ш = ш(я) в круге 14 = (г! я! < 1) однолистна при а! < 1/4, и поэтому она же отображает круг Ь на область ??.
'1игленно конформное отображение строилось для значений параметра а = 0,1! 0,2: 0,24. Заметим, что а = 0,25 — критическое значение параметра и что контур, ограничивающий областгь имеет угловые точки. Обратим вник!ание читателя на характер убывания величины е в первых двух колонках. Такое поведение погрешности типично для ненасыщаемых методов.
Е'ис. 1. Кривая и! = (с+0,2я~! .=-ехр(10), 0~ (0 ~ (2я? 419 з 5, Кубаглурвыс формулы Рис. 2. Кривая ю = 1я + 0,24 л: з =- ехр11В1, 0 < В ( 2я1 На рис. 1, 2 изображена область т1, отвечаюшая значениям параметра а = 0,2; 0,24. Хотя прн а =. 0,24 мы не имеем угловых точек па контуре, но максимум кривизны равен 1775, Таким образом, формально граница гладкая и даже аналитическая, но производные от функции г(В) огромные, и для того, чтобы ликвидировать зту неблагоприятную особенность, нужны соответствующие значения параметра п. Так, возможно, что расчет при а = 128 —: 140 вполне удовлетворит по точности самые взыскательные требования.
Глйвл 7 Численное решение задачи Коши для обыкновенных дифференциальных уравнений Мы рассмотрели вопросы приближенного вычисления определенных интегралов; полученные при этом результаты могут быть использованы полностью при вычислении неопределенных интегралов или же при численном решении обыкновенных дифференциальных уравнений. Многие формулы численного решения обыкновенных дифференциальных уравнений получаются на основе тех или иных квадратурных формул.
Однако в задаче о неопределенных интегралах или в решении задачи Коши имеется существенно новый элемент, который обусловливает отличие этих задач от задач о вычислении определенных интегралов. В самом деле, при вычислении неопределенного интеграла или решении дифференциального уравнения мы имеем возможность использовать ранее найденные значения искомого решения. Это нам позволяет либо нолучить формулу более высокого порядка точности, либо сделать сам алгоритм проще.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
