Федоренко Р.П. Введение в вычислительную физику (Федоренко Р.П. Введение в вычислительную физику.djvu), страница 5

Но величина не- вязки г = [(г'г'х1)~+ ( р/хз) [н~ и, следовательно, шаг з такой инвариантностью не обладают. При этом возникает проблема выбора «праенльных» масштабов х„хи В своей практике автор в подобных ситуациях руководствовался правилом, условно названным «принципом равноправия»: масштабы нужно выбирать такими, чтобы одинаковые изменения х н у приводила к численно близким изменениям У и р. Формулы для малых приращений у и р показывают, что эта цель в известной мере будет достигнута (в окрестности данной точки (х, у)), если взять х = (г'з + г»)11» х = (~рз + <рз)11з Таким образом, мы приходим к модифицированному методу Ньютона с нормировкой.

Алгоритм стандартного шага в точке (х, у) дополняется следующим: после вычисления производных и направления (Ьх, Ьу) вычисляются «масштабы» хп и шаг з выбирается минимизацией масштабированной невязки. Эффект этого прн- ОСНОВЫ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ 1ч. ! ема иллюстрирует табл, 2. Обозначения в ней — те же, что и в табл. 1, только величины гв и г, означают величины масштабнрованной невязки в точке (хь, уа) и в следующей точке (ха+!, ух +!). Заметим, что теперь у нас нет единой невязки, которая убывала бы в процессе решения, чем в сущности и обеспечивается сходимость метода. В точке (х«+!, у«+!) есть две невязки! при масштабировании в точке (хь, уь) и при масштабировании в точке (х"+', уь+!). Теорема о сходимости модифицированного метода Ньютона утрачивает силу, но зато сама сходнмость стала существенно лучше.

Данные г" =0 и р 0 в последней строке табл. 2 означают, что зти величины не больше 5.10 6. На основании вышесказанного естественно возникает вопрос: действительно ли зто малая величина? Ответ на него несложен. Если вычислитель постулировал„что для Таблица 2 го 5.72 12.75 2,00 0.60 0.33! 0.03! О.оооб 5.52 1 1.59 0.55 0.0826 0.026 0.00054 0 0.066 О.! 25 отш 1.90 2.00 1.0 1.0 х и у величина 10 5 является малой, то естественно считать малыми для г и !р изменения, порождаемые такими малыми Ьх и Ьу. Это опять-таки приводит нас к тому масштабированию, которое было использовано, и величины г, р порядка 5.10 6 следует тоже считать малыми.

Из табл. 2 видно, что в единицах к„к величина / = 111 стала «малой» по сравнению с р = 17. Допускается ее существенное увеличение ради уменьшения !р. Внешне большое значение У = 338 затем сравнительно малыми изменениями Ьх и Ьу доводится до нуля. Это есть следствие разной чувствительности 5 и р к изменениям х, у, т.е. существенно разных величин их производных.

После того как два-три раза подряд минимизация г по В в модифицированном методе Ньютона приводит к значениям, близким к единице, переходят на обычный метод Ньютона, не тратя машинного времени на подбор 5. Однако в данном примере, после получения 2.0000 2.414 3.209 2.89! 1.806 1.157 !. ! 4265 1.14220 3.0000 2.626 0.542 -0.389 -0.578 -0.460 -0.48663 -0.48626 ! ! !.Ооо 127.500 338.4 200,0 17.3 0.1! 5 0.004 о 17.000 16.32 ! 4,66 3.07 -1.! 4 0.516 -0.006 о к!1 гкшкннк систем нклннкйных м чзнкннй приближения, хорошего в смысле теоремы 1, сходимость оказывается столь быстрой, что этот прием не дал бы нам никакой экономии. Другая, но по существу близкая, нормировка была предложена немецким математиком Ддфлхардом.

В его варианте модифицированного метода Ньютона после вычисления г(х) и матрицы Р,(х) вычисляется матрица А =г,' и минимизируется невязка г~(к) = (А У(х+ к Ьх), А У(х+ к Ьх)). Смысл этой конструкции станет ясен, если проанализировать поведение правой части при малых Ь (в первом приближении): АГ(х+ Ь) = АГ(х) + А Г„(х) Ь= А У(х) + Ь.

Таким образом, в окрестности точки х невязка устроена очень просто: '(Ь) =С+(в, Ь)+(Ь, Ь), Такая ситуация наиболее благоприятна для алгоритмов поиска минимума, а решение системы г'(х) = 0 можно трактовать как поиск минимума гк ш ~~у(х) ~)к, Метод продолжения по параметру. Опишем в общих чертах другой прием, имеющий ту же цель — ослабить требования к выбору начального приближения и обеспечить надежную сходимость метода решения системы уравнений г (х) = О. В литературе этот метод иногда называют «методом инвариантного погружения». Рассмотрим семейство задач Р(х, Х) = О, где Х вЂ” скалярный параметр.

Сконструируем это семейство так, что Р(х, 1) к«У(х), а при А = 0 уравнение Р(х, 0) = 0 легко решается или даже имеет явное решение. Это и есть «погружение» исходной задачи в семейство задач. Формально построить такое погружение часто не представляет труда. Вот, например, самый простой способ: Р(х, Х) ш (1 — Х)х+ Хг(х). (5) Уравнение Р(х, 0) = 0 имеет очевидное решение х =О. Пусть х( Х) есть решение уравнения (5). Последовательно решается серия задач.

Имея решение х(Х), меняем А на Х+ йА и решаем уравнение Р(х, Х + ЬА) = О тем или иным итерационным методом. Используем х(А) как хорошее начальное приближение (ЬХ, естественно, считаем малым изменением). Таким образом можно (прн благоприятном ходе событий) добраться до Х = 1 и получить решение исходной задачи. Нетрудно оформить эту конструкцию в виде задачи Коши для системы обыкновенных дифференциальных уравнений, т.е.

вычислить производную и'хИА. В самом деле, дифференцируя по А, получаем 0= — Р(х, й) =Р„(х, Х) — "„+ Р„(х, А), (ч.( ОСНОВЫ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ откуда †" = — Р„((х, Х) Р (х, Х). Связь с методом Ньютона достаточно прозрачная. Сведение к задаче Коши иногда считают исчерпывающим решением проблемы, поскольку многие полагают эту задачу самой простой в вычислительной математике. Это мнение (ошибочное в столь общей форме, как мы увидим в дальнейшем) основано на том факте, что для решения задачи Коши существуют не только строго обоснованные алгоритмы, но даже стандартные программы н можно, не имея представления о том, как они работают, просто обращаться к ним.

Реализуя метод продолжения по параметру, часто сталкиваются с тем, что график х(Х) имеет 5-образную форму, т.е. имеется несколько ветвей решения уравнения (г(х, А) = О. Отслеживам одну из них, достигают некоторой точки А(, за которую даннам ветвь х(А) не продолжается, — решения уравнения Г(х, А+ (ьА) = О, близкого к х(А), не существует. Внешне это проявляется в вырождении матрицы г"„, т.е. де( Г„-» О. Для продолжения решения задачи нужно двигаться по Х в обратном направлении, перейдя, однако, на другую ветвь функции х(А). Технически зто реализуется следующим образом: л-мерный вектор х и скалярный параметр Х обьединяются в единый (л + 1)-мерный вектор у = (у„..., у„,).

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

Признаком приближения точки поворота длм компоненты у,. может служить, например, уменьшение ее приращения за один шаг. В 2, Численное дифференцирование Практическое применение приведенных в й 1 формул численного дифференцирования связано с необходимостью выбора подходящего шага л. Возникающие здесь проблемы рассмотрим длм простоты на примере функции талька одного переменного У(х). Нас интересует погрешность формулы численного дифференцирования. Тривиальный («школьный») ответ «чем меньше л, тем точнее формула численного дифференцирования» основан на известном соотношении ~Ч ) У(х+Ю-У(х) численное днфрегенцигоелние Однако, как мы увидим, при реальных вычислениях ситуация сложнее.

Пусть функция /(к) гладкая, но, работая на ЭВМ, мы имеем дело не с /(х), а с ее машинным представлением /"(х) = /(х)(1 + а(х)), где а(х) — относительная погрешность вычисления /. Разумеется а зависит от к, но для всех интересующих вас х пусть имеется оценка ~ а(х) ~ в а м 1. Величина а может быть связана хотя бы с конечным числом знаков в представлении / в памяти Э ВМ (аж10 )а на БЭСМ-б, а 10 ~ на ЕС, ащ10 и на ЕС при двойной точности). Но если /(х) вычисляется достаточно сложно, погрешнос)ь а может быть и существенно большей величиной, не всегда допускающей хорошую оценку. Таким образом, используя численное дифференцирование, вычисляем /"(х+Л)-/ са) /(х+Л) — /(х) + а~ а) Л Л / Л Разложим в ряд Тейлора: /(х + Ь) = /+ Ь /„+ — "~~ /„„+ О(Ьз), Итак / (х+Л) — / (х) / ( ) 1 Л / 1 О(Ьа) + О ы / Погрешность численного дифференцирования состоит из двух частей.

Первая из них связана с заменой оператора дифференцирования оператором конечной разности. Она имеет величину О(/„„Ь) и стремится к нулю при Ь - О. Степень Ь в погрешности аппроксемации называют порядком аппроксимации Вторая часть погрешности связана с неточностью вычисления /. Она имеет величину О(а//Ь) и при Ь- 0 стремится к бесконечности. Полная погрешность численного дифференцирования есть сумма погрешностей аппроксимации и округления: 0.5Ь/„„+ 2а//Ь. Легко вычислить шаг Ьр, при котором полная погрешность мини.о~ д,а, п)еп,,г.э рн ~...: ~, ж, числяется через а и /„„, точные значения которых обычно неизвестны.

Существуют достаточно надежные алгоритмы численного дифференцирования, основанные на вычислении большего числа значений /. По этим значениям вычисляются варианты разноспюй производной, из сравнений которых между собой отбирается наиболее достоверное значение. Надежность подобных алгоритмов оплачивается большим объемом вычислений. Часто, например при реше- основы еычнслнтгльноя млтемлтнкн [Ч. ( 26 нин уравнений методом Ньютона, нет нужды в особенно высокой точности численного дифференцирования н не требуется очень точно определять 6 .

В частности, автор иногда использовал простой способ проверки того, является ли данное Ь подходящим для численного дифференцирования. Этот способ основан на подсчете числа «сокращающихся знаков». Поясним его иа примере численного дифференцирования функции сх в точке х =!. Пусть сх вычисляется с Таблица 3 семью верными знаками. Вычислим ее разностную производную с шагом Ь, равным 0.1, 0.01, 0.00!, 0.0001. В табл. 3 пояснения требует только последняя строка: Л вЂ” это число сократившихся главных знаков, Таким образом, наилучший результат получается при сокращении половины знаков. В общем случае число «сократившихся знаков» при численном дифференцировании можно оценивать величиной ) У(хч-Ы -У(х) ! е гехтет и ~ ~ Заметим, что в этом простом примере мы сталкиваемся с одной из самых грозных опасностей в приближенных вычислениях: если результат получается при вычитании двух очень близких друг к другу величин, происходит резкая потеря точности, относительная погрешность результата сильно возрастает.

Это явление носит название «сокращение знаков» и доставляет массу неприятностей. Оценим погрешность формулы центральных разностей: «у у" (х -ь Л) — у" (х — Л) у(х+ Ы вЂ” у(х — Л) ) е ах 2Л 2Л + О ~Л ~) . Разлагая в ряд Тейлора У(х + Ь) и /(х — Л), получаем 2 (х+Л)-/ (х — Л) у т )+Л у +Отав)+О зг) численное диоовгвннигооАниз Главный член погрешности аппроксимации имеет второй по /( поря- док, а оптимальный шаг Таким образом, формула центральной разности точнее формулы односторонней разности, но требует большего шага й (при одинаковых шагах сокращается больше знаков). Например, еьоо( еоооо 2 721001 2 715565 = 0 005436, (ех)„ж 2.718 (т.е.

получается та же точность, что и в формуле односторонней разности). При /з = 0.01 имеем е'о' — еото = 2.745601 — 2.691234 = 0.054367, (е")„2.71835. В дальнейшем нам часто придется использовать формулу приближенного вычисления второй производной Ез/ /(х+А) — 3/(х)+/(х-О) о хо дз Оценим точность этой формулы, Очевидно, погрешность округления есть (4о//(з)/. Вычислим погрешность аппроксимации, используя разложение в ряд Тейлора /(х + Л) н /(х — /)). После простых вычислений найдем полную погрешность численного дифференцирования: дз н оптимальный шаг: /з () 48 о///(~ () п~ Вообще, из всех полученных выше формул для Ао видно, что «наилучший шаг» тот, при котором погрешности аппроксимации н округления совпадают (близки).