УМК (1013374), страница 23
Текст из файла (страница 23)
Если hi +1 = h = const , сетка называется равномерной (регулярной), а если hi +1 = var – неравномерной (нерегулярной). В случае равномер-ной сетки узлы находятся по формуле x i = x 0 + ih, i = 0, n .Решение находится в виде последовательности значений yˆ0 , yˆ1 , yˆ2 ,..., yˆn , являющихся приближением значений y 0 , y ( x1 ), y ( x 2 ),..., y ( x n ) точного решения y(x ) в узлахсетки Ω n (рис. 1).135yŷ iy ( xi −1 )y ( x1 )y ( xi +1 )ŷ i −1y ( xi )y = y(x )ŷ i +1ŷ1ŷ 0y ( xn )y00ŷ nx0x1xi −1xixi +1xnxРис. 1Численные дискретные методы решения обыкновенных дифференциальных уравнений, позволяющие найти решение только в узлах сетки, делятся на две группы: явные инеявные.Значение ŷi +1 на (i + 1) -м шаге может определяться явно:yˆi +1 = Φ( x i − k +1 ,..., x i −1 , x i , yˆi − k +1 ,..., yˆi −1 , yˆi ) ,где Φ(.) – некоторая функция, зависящая от конкретного метода (кроме последней рассчитанной точки ( x i , yˆi ) могут использоваться еще (k − 1) предыдущих точек), илинеявно:yˆi +1 = Φ( x i − k +1 ,..., x i −1 , x i , x i +1 , yˆi − k +1 ,..., yˆi −1 , yˆi , yˆi +1 ) ,где искомая величина ŷi +1 входит одновременно и в левую, и в правую часть.Явные и неявные методы делятся также на одношаговые и многошаговые ( k шаговые).
В одношаговых методах для расчета очередной точки ( x i +1 , yˆi +1 ) требуетсяинформация только о последней рассчитанной точке ( x i , yˆi ) . В k -шаговых методах длянахождения точки ( x i +1 , yˆi +1 ) требуется информация о k предыдущих точках.Формулы явных или неявных методов в общем случае представляют собой нелинейные уравнения относительно ŷi +1 и называются разностными схемами.Локальной ошибкой численного метода на (i + 1) -м шаге называется величинаε i +1 (h) = yˆi +1 − y ( x i +1 ) ,где y ( x i +1 ) – значение точного решения при x = x i +1 , а ŷi +1 – приближенное решение,получаемое по формулам при условии, что вместо приближенных значенийyˆi , yˆi −1 ,..., yˆi − k +1 используются значения, соответствующие точному решению, т.е.y ( x i ), y ( x i −1 ),..., y ( x i − k +1 ) .136Глобальной ошибкой называется величина e n (h) = yˆn − y ( x n ) , где ŷ n – значение,получаемое по формулам при i = n − 1 .Глобальная ошибка определяется:а) ошибками округления и ошибками арифметических действий, обусловленнымичислом разрядов компьютера и характером выполняемых операций для расчета значенияискомой функции в очередной точке x i +1 ;б) методическими ошибками, определяемыми выбранным алгоритмом;в) переходными ошибками, обусловленными тем, что при расчете значения ŷi +1вместо точных значений y ( x i ), y ( x i −1 ),..., y ( x i − k +1 ) берутся приближенные значенияyˆi , yˆi −1 ,..., yˆi − k +1 , полученные на предыдущих шагах.Локальные ошибки «переносятся» в точку x n и формируют глобальную ошибку.Число p называется порядком (точностью) численного метода, если его глобальная ошибка есть О большое от h p , т.е.
e n (h) = O (h p ) .Пояснение. Пусть R (h) – некоторая функция переменной h (как правило, R (h) –остаточное слагаемое некоторой аппроксимационной формулы) с конечной областью определения DR на полуоси h > 0 , причем h ∈ DR . Тогда, если при некотором h ≤ h0справедливо неравенство R (h) ≤ ch k , где c = const , не зависящая от h , k – целое число, h0 > 0 , то пишут R (h) = O (h k ) и говорят, что R (h) есть «O большое от h k » приh → 0.На практике в качестве характеристики точности метода часто используется велиyˆi − y ( x i ) .чина ε(h) = maxi = 0,1,..., nМожно показать, что если локальная ошибка имеет порядок ( p + 1) , т.е.ε i +1 (h) = O (h p +1 ) , то глобальная погрешность имеет на единицу меньший порядок, т.е.e n (h) = O (h p ) .Перейдем теперь к рассмотрению устойчивости численных методов.
Она проверяется на «тестовом примере»y ′ = μ y,y (0) = 1 ,где μ – в общем случае комплексная константа. Дифференциальное уравнение являетсяпростейшим линейным уравнением, и для него можно получить значимые критерии устойчивости в явной форме.Метод называется ограниченно устойчивым, если существует такое число hкр. > 0 ,что при использовании метода для решения тестового примера, где Re μ < 0 , с шагом0 < h < hкр.
при i → ∞ глобальная ошибка ограничена. Величина hкр. называется критическим шагом. Если h > hкр. , глобальная ошибка может неограниченно возрастать. В ограниченно устойчивых методах при задании величины шага h необходимо учитыватьзначение критического шага hкр. . Для сложных дифференциальных уравнений и систем137нахождение hкр. является самостоятельной задачей, а свойство ограниченной устойчивости предупреждает вычислителя о возможных проблемах. Поэтому на практике становится актуальной задача конструирования таких методов, которые были бы устойчивыпри любом значении шага, а его величина выбиралась бы только исходя из желаемойточности расчетов (при этом класс решаемых задач может быть ограничен).Метод называется А-устойчивым, если при его применении с любым фиксированным положительным шагом h все численные решения тестового примера с комплекснойконстантой μ ( Re μ < 0 ) стремятся к нулю при i → ∞ .А.
ЯВНЫЕ МЕТОДЫ РЕШЕНИЯ СИСТЕМ ОБЫКНОВЕННЫХДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙА1. Явный метод ЭйлераРассмотрим проблему нахождения численного решения задачи Коши:dy= f ( x , y ),dxy(x 0 ) = y0 .Вводится в общем случае неравномерная сетка Ωn = (x0 , x1,...,xi , xi +1,...,xn ) . Величина шага hi +1 = x i +1 − x i выражается через узловые точки. Для аппроксимации производной⎛ dy ⎞⎜ ⎟⎝ dx ⎠x = xiиспользуем⎛ dy ⎞Ш 2,i = (x i , x i +1 ) : ⎜ ⎟⎝ dx ⎠x = xi=формулу,y i +1 − y ihi +1записанную+ O (hi +1 )надвухточечномшаблоне⎛ h i +1⎞M 2,i ⎟⎟ .⎜⎜⎝ 2⎠Далее заменяется правая часть уравнения ее сеточным представлением, т.е.f ( x , y ) → f ( x i , yi ) , а вместо y ( x ) рассматривается сеточная функция yˆi ≈ y ( x i ) , которая определяется только в точках сетки.
Выполняется подстановка аппроксимаций производной и правой части в дифференциальное уравнение:y i +1 − y ihi +1+ O(hi +1 )= f ( xi , y i ) .После отбрасывания остаточных слагаемых получается явная схема Эйлера первого порядка (явный метод Эйлера):yˆi +1 = yˆi + hi +1 f ( x i , yˆi ) ,i = 0, n − 1 , yˆ0 = y 0 ;Порядок точности метода, как правило, определяется порядком аппроксимациисхем, явный метод Эйлера является ограниченно устойчивым с критическим шагом2hкр.
= − (см. тестовый пример).μ138А2. Метод Эйлера-КошиДля аппроксимации производной применяется формула:⎛ dy ⎞⎜ ⎟⎝ dx ⎠x = xi=⎛ h2⎞⎜M 3,i ⎟ .⎜ 6⎟⎝⎠y i + 1 − y i −1+ O (h 2 )2hВыполняется подстановка аппроксимаций производной и правой части в дифференциальное уравнение:y i +1 − y i −12h+ O(h 2 ) = f ( x i , y i ) .После отбрасывания остаточных слагаемых получается явная схема метода Эйлера–Коши второго порядка:yˆi +1 = yˆi −1 + 2h ⋅ f ( x i , yˆi ) , i = 1, n − 1 .Для начала расчетов требуется иметь две «разгонные» точки yˆ0 , yˆ1 . Первая определяется известным начальным условием yˆ0 = y 0 , а вторая может быть найдена с помощью другого метода, например, по формуле: yˆ1 = y 0 + h1 f ( x 0 , y 0 ) .А3.
Модифицированный метод ЭйлераМодифицированный метод Эйлера второго порядка:yˆi+12= yˆi +hi +12f ( x i , yˆi ) ,i = 0, n − 1 ,h⎛⎞yˆi +1 = yˆi + hi +1 f ⎜⎜ x i + i +1 , yˆi + 1 ⎟⎟ ,22 ⎠⎝i = 0, n − 1 .Интервал устойчивости h μ ∈ (−2, 0) (здесь μ – действительное число в тестовомпримере) модифицированного метода Эйлера совпадает с интервалом устойчивости явного метода Эйлера.А4. Методы Рунге-КуттыФормулы семейства методов Рунге–Кутты имеют следующую структуру:[]yˆi +1 = yˆi + hi +1 ⋅ b1K 1,i + b2 K 2,i + " + bs K s ,i ,139yˆ0 = y 0 ,i = 0, n − 1 ,K 1,i = f ( x i , yˆi ) ;гдеK 2,i = f ( x i + c 2 hi +1 , yˆi + hi +1 a2,1K 1,i ) ;K 3,i = f ( x i + c3 hi +1 , yˆi + hi +1 (a3,1K 1,i + a3,2 K 2,i ) ;#K s ,i = f ( xi + c s hi +1 , yˆi + hi +1 (as ,1K 1,i + as ,2 K 2,i + " + as , s −1K s −1,i ) ,где s – число стадий (этапов), K s ,i – значения коэффициентов схемы Рунге–Кутты, вычисленныенаосновеправойчастидифференциальногоуравнения,c j , j = 2, s; al , m , l = 2, s ; m = 1, s − 1; bk , k = 1, s .
Первый индекс в обозначениях коэффициентов является порядковым номером, а второй соответствует индексу точки x i –началу отрезка [ x i , x i +1 ] , на котором производится расчет.В некоторых методах кроме вычисления приближенного решения ŷi +1 определяyi +1 по формулеется еще дополнительное значение ~[~yi +1 = yˆi + hi +1 ⋅ b~1K 1,i + b~2 K 2,i + " + b~s K s ,i],порядок которого, как правило, на единицу больше или меньше обеспечиваемого выраyi +1 служит для учета погрешности и управления вежением для ŷi +1 . Величина yˆi +1 − ~личиной шага.Наибольшее распространение в вычислительной практике нашел метод Рунге–Кутты четвертого порядка:yˆi +1 = yˆi +гдеhi +16(K 1,i + 2K 2,i + 2K 3,i + K 4,i ),yˆ0 = y 0 , i = 0, n − 1 ,hh⎛⎞K 2,i = f ⎜⎜ x i + i +1 , yˆi + i +1 K 1,i ⎟⎟ ,22⎝⎠hh⎛⎞K 3,i = f ⎜⎜ x i + i +1 , yˆi + i +1 K 2,i ⎟⎟ , K 4,i = f x i + hi +1 , yˆi + hi +1 ⋅ K 3,i .22⎝⎠Схема является четырехчленной, первый коэффициент K 1,i относится к точке x i ,K 1,i = f i = f ( x i , yˆi ),(второй и третий – к средней точке x i +выбираютсяследующиеhi +12), четвертый – к точке x i +1 .
Для этой схемыпараметры:s = 4,c 2 = c3 =1; c4 = 1 ;211; a3,1 = a4,1 = a4,2 = 0 ; a3,2 = ; a4,3 = 1 .22Метод Рунге–Кутты, как и методы Эйлера, является одношаговым, так как значение yˆ i +1 вычисляется на основе текущего значения yˆ i . По сравнению с явным методомЭйлера здесь на одной итерации требуется вычислять значение правой части решаемогоуравнения четыре раза. Как и явный метод Эйлера, метод Рунге–Кутты не требует дополнительных разгонных точек, что позволяет легко менять шаг в процессе вычислений.В методе Рунге–Кутты пятого порядка точности для расчета точки yˆ i +1 используются следующие соотношения:a2,1 =140yˆ i +1 = yˆ i +гдеhi +16( K 1,i + 4K 3,i + K 6,i ) ,yˆ0 = y 0 , i = 0, n − 1 ,hh⎛⎞K 2,i = f ⎜ x i + i +1 , yˆ i + i +1 K 1,i ⎟ ,22⎝⎠K 1,i = f ( x i , yˆ i ),hh⎛⎞K 3,i = f ⎜ x i + i +1 , yˆ i + i +1 ( K 1,i + K 2,i ) ⎟ , K 4,i = f ( x i + hi +1 , yˆ i − hi +1 K 2,i + 2hi +1 ⋅ K 3,i ) ;24⎝⎠2hh⎛⎞K 5,i = f ⎜ x i + i +1 , yˆ i + i +1 (7 K 1,i + 10 K 2,i + K 4,i ) ⎟ ,327⎝⎠hh⎛⎞K 6,i = f ⎜ x i + i +1 , yˆ i + i +1 (28K 1,i − 15 K 2,i + 546 K 3,i + 54 K 4,i − 378 K 5,i ⎟ .25625⎝⎠А5.
Методы Адамса–БэшфортаМногошаговые схемы Адамса–Бэшфорта:– второго порядка:hyˆi +1 = yˆi + [3 f i − f i −1 ] , i = 1, n − 1 ;2– третьего порядка:hyˆi +1 = yˆi + [23 f i − 16 f i −1 + 5 f i − 2 ] , i = 2, n − 1 ;12– четвертого порядка:hyˆi +1 = yˆi +[55 f i − 59 f i −1 + 37 f i − 2 − 9 f i − 3 ] , i = 3, n − 1 ;24– пятого порядка:yˆ i +1 = yˆ i +h[1901 f i − 2774 f i −1 + 2616 f i −2 − 1274 f i −3 + 251 f i −4 ] , i = 4, n − 1 ,720где f i = f ( x i , yˆ i ) .Для начала расчетов по первой формуле требуются две «разгонные» точки: yˆ0 , yˆ1 ,по второй формуле – три «разгонные» точки: yˆ0 , yˆ1 , yˆ2 , по третьей формуле – четыре«разгонные» точки: yˆ0 , yˆ1 , yˆ2 , yˆ3 , по четвертой формуле – пять «разгонных» точек:yˆ 0 , yˆ 1 , yˆ 2 , yˆ 3 , yˆ 4 . Их необходимо вычислить с порядком точности не меньше порядка точности схемы.Методы Адамса–Бэшфорта не позволяют изменять шаг в процессе расчетов.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
