1626435587-51311eae4652e8ad616b5bdef025cbb3 (844239), страница 20
Текст из файла (страница 20)
Если неустойчивость проявляется при использовании настолькобольших шагов сетки ℎ, что получаемое численное решение уже не аппроксимирует точное решение ОДУ, зачем в таком случае нужно вводить понятие условной устойчивости? Зачем вообще рассматривать иисследовать численные решения на сетках со столь большим шагом ℎ?Ответ становится понятным при переходе от одного уравнения (70)к системе ОДУ. Рассмотрим суть проблемы на самом простом примере — системе двух независимых линейных ОДУ:{︂ ′{︂1 = −11 () = 1 (0) · −⇒(94)′2 = −100022 () = 2 (0) · −1000Функция 2 () быстро затухает, так что поведение решения u() =(1 (), 2 ()) практически полностью определяется первой компонентой 1 ().
Тем не менее, шаг численного интегрирования системы ОДУ102должен определяться компонентой 2 , несущественной с точки зренияфизики. Для обеспечения устойчивости численного решения u(), полученного методом Эйлера, в соответствии с (92) необходимо выбратьℎ ≤ 0,002, в противном случае 2 будет экспоненциально возрастать —возникнет неустойчивость. В следующих лекциях будет показано, чтоданная проблема усугубляется при численном решении уравнений вчастных производных, разностные схемы для которых являются, посути, системой большого количества ОДУ.Система (94) называется, поскольку имеет существенноразличные (в частном случае (94) отличающиеся в тысячу раз) масштабы изменения различных компонент решения 1 и 2 (в общемслучае — некоторых их линейных комбинаций).
Обобщим сказанное наслучай нераспавшихся систем.. Говорят, что система линейных уравнений u′ = f, если велико отношение максимального и минимального модуля собственных значенийжёсткойОпределениеляется жёсткой=яв-max Re| |,min Re| | < 0, = 1, 2, . . . , .числом жёсткости системыУказанное отношение называют. В случае, если матричные коэффициенты зависят от , система может оказаться жёсткой в некоторой области по .В качестве примера рассмотрим систему из двух уравнений с постоянными коэффициентами [4],[A11]:{︂ ′998 + 1998 , =(95) ′ = −999 − 1999 .Собственные числа матрицы системы (95) равны 1 = −1, 2 = −1000,откуда число жёсткости равно 2 /1 = 103 ≫ 1 и, следовательно, система (95) является жёсткой.
Общее решение системы имеет вид(︂)︂(︂)︂(︂)︂()21=·− + ·−1000 .()−1−1При ≈ 1, ≪ значения в правой части (95) будут порядка 1, однакопри использовании метода Эйлера (73) следует выбирать шаг сеткиℎ < 0,002, в противном случае численное решение будет неустойчивым.Использование неявных схем для численного решения систем ОДУпозволяет существенно повысить скорость вычислений за счёт использования сетки с меньшим числом узлов. Очевидно, что за это приходит103ся платить необходимостью решения конечного (чаще всего, нелинейного трансцендентного) уравнения на каждом шаге численного интегрирования.
Рассмотрим данную задачу на примере простейшей неявной схемы первого порядка (93). Для выполнения шага численного интегрирования нам необходимо разрешить уравнение (93) относительно+1 . Какими методами может быть решена эта задача?Метод дихотомии не обобщается на системы уравнений и потомуне подходит для данной задачи. Посмотрим, возможно ли применениеметода простых итераций и метода Ньютона — Рафсона.Для использования итерационных методов необходимо, чтобы итерационный процесс сходился, и, кроме того, нужно указать условиезавершения итераций при вычислении +1 в программе. На практике зачастую пользуются так называемымиметодами,которые заключаются в выполнении фиксированного числа итераций(например, одной или двух). Несложно понять, что при использовании конечного фиксированного числа шагов результат использованияитерационного процесса — значение +1 — можновыразить через .
Как следствие, полунеявные методы по сути являются разновидностью явных численных схем и потому не могут обладатьабсолютной устойчивостью. Тем не менее, они успешно используютсяпри интегрировании жёстких систем, позволяя существенно ослабитьограничение на шаг сетки и достичь тем самым выигрыша в быстродействии.Рассмотрим использование неявных схем на примере системы издвух уравнений общего вида:{︂ ′ = (, ),(96) ′ = (, ).полунеявнымиявным образомЗависимость и от мы предполагаем, но для краткости не пишем.Будем решать систему (96) по неявной схеме (93). Полагая функции и гладкими, можем записать:{︂ ≡ +1 − = ℎ + ℎ · + ℎ · + (2 , 2 ), ≡ +1 − = ℎ + ℎ · + ℎ · + (2 , 2 ),где для краткости значения функций и их частных производных на( , ), . .
. Удобно-м шаге по обозначены ≡ ( , ), ≡ переписать полученные соотношения в матричном виде:(︂)︂(︂)︂= ℎˆ+ ℎf + ( 2 , 2 ),(97)(︂)︂(︂)︂ ˆ ≡, f≡. 104Ввиду наличия нелинейных слагаемых ( 2 , 2 ) уравнение (97) в общем случае не решается точно. Попытаемся построить его численноерешение, используя вначале метод простых итераций:)︃)︃(︃(︃()(+1)+ ℎf + ( 2 , 2 ).(98)= ℎˆ()(+1)здесь верхний индекс обозначает номер итерации, нижний — узел сетки. Исследуем сходимость итерационного процесса (98). Пусть u* — решение системы (97), т.
е. в первом порядке по ℎ можно записатьˆ * + ℎf .u* = ℎuДля того, чтобы итерационный процесс сходился, необходима устойчивость по малым отклонениям от u* . Подставляя в (97) ( () , () ) =u* + () , получаем с точностью до членов второго порядка малостиˆ () . (+1) = ℎСледовательно, итерационный процесс (98) будет сходиться, если всесобственные числа матрицы ℎˆ (97) по модулю меньше 1.
(В случае одной переменной данное условие переходит в |′ | < 1 — условие сходимости итерационного процесса (+1) = (() ) для уравнения = ().)Указанное условие приводит нас к ограничению на максимальное значение ℎ, аналогичное условию (92). Действительно, как и (92), условиеℎ · | ()| < 1 накладывает ограничение на произведение шага сетки ℎна некоторую комбинацию частных производных , , и , входящих в матрицу Якоби ˆ.Заметим, что уменьшение ℎ при поиске +1 хотя и обеспечиваетсходимость итерационного процесса (98), но при этом приводит ровно к той же потере скорости счёта, что и ограничение на шаг сетки,связанное с условием устойчивости (92) для явной схемы Эйлера (73).В этой связи представляется целесообразным использовать для численного решения уравнений (97) метод Ньютона.
По сути, шаг итераций в методе Ньютона есть точное решение линеаризованной задачи(97). Решая (97) в линейном приближении по ℎ, имеем:[︁]︁ (︂ )︂ˆ − ℎ= ℎf .Организуя итерационный процесс в соответствии с полученными формулами и записывая номер итерации в верхнем индексе:(︃)︃(+1)[︁]︁−1 (︁)︁()= − ℎˆ ()ℎf () + u − u+1 ,(+1)105(︃(+1)+1(+1)+1(︃)︃=()+1()+1)︃[︁+ − ℎˆ()]︁−1(︃() − +1 + ℎ ()() − +1 + ℎ ())︃(99)где () , () и ˆ() — значения функций и их частных производных,()()вычисленные в точке (+1 , +1 ).
В качестве начального приближения(0)(0)(+1 , +1 ) в (99) естественно использовать ( , ).В случае линейных систем с постоянными коэффициентами (таких,например, как (94) и (95)), однократное применение формулы Ньютона(99) даёт точное решение задачи: итерационный процесс (99) сходитсяза один шаг.В общем случае cходимость итерационного процесса (99) не гарантируется, однако в ряде случаев данный метод может быть успешноиспользован на практике, обеспечивая приемлемую точность при использовании 1-2 итераций для неявной схемы первого порядка (93).Как и в случае явных схем, в расчётах целесообразно использоватьметоды более высокого порядка точности при условии, что входящиев уравнения функции являются достаточно гладкими.
Неявные схемы, возникающие в результате обобщения методов Рунге — Кутты, известны как методы Розенброка (Rosenbrock), а также методы Канса —Рентропа (Kans, Rentrop). Неявные многошаговые методы Адамса известны как методы Адамса — Мултона (Adams, Moulton).6.8. ЗаключениеТаким образом, в данной главе был рассмотрен ряд базовых методов интегрирования обыкновенных дифференциальных уравнений и их106систем.
Наиболее простым для анализа и реализации является явныйметод Эйлера, однако на практике следует избегать его использованияввиду крайне низкой скорости сходимости. Наиболее употребительнымв приложениях является метод Рунге — Кутты 4-го порядка точности.Достичь ещё более высокого порядка точности (что имеет смысл лишьдля достаточно гладких функций в правой части ОДУ, которые притом хорошо интерполируются полиномами) позволяют методы Адамса. К их недостаткам относится отсутствие самостарта, относительнобольшие численные коэффициенты в остаточном члене и отсутствиеабсолютной устойчивости даже при использовании неявных модификаций многошаговых методов (см. [A8, с. 357] и [A1, c. 277]). При интегрировании жёстких систем ОДУ численные решения, построенные появным схемам, могут быть неустойчивыми при недостаточно мелкомшаге сетки.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
