Блюмин А.Г., Федотов А.А., Храпов П.В. - Численные методы (1078410), страница 6
Текст из файла (страница 6)
Задача Коши для обыкновенных дифференциальных уравненийОпределение 3. Численный метод решения задачи (3.1.1) называется методом k-ого порядка точности, еслиε(m) = O(hk ) при h → 0,или иначеε(m) 6 C · hk при h → 0,(3.2.5)где C = const > 0 зависит от f (x, u), a, b и численного метода решениязадачи (3.1.1), но не зависит от h.Покажем, что метод Эйлера является методом первого порядка точности. Будем полагать, что арифметические вычисления проводятся точно ипоэтому локальная и методическая ошибки совпадают. Докажем следующееутверждение.Лемма. Пусть при любом i = 0, 1, ..., n − 1 справедлива оценкаεi+1 6 (1 + c1 · h)εi + c2 · hk+1 ,(3.2.6)где h = (b − a)/n; c1 > 0 и c2 > 0 — постоянные, не зависящие от h. Тогдасоответствующий метод численного интегрирования задачи (3.1.1) являетсяметодом k-ого порядка точности.Доказательство.
В соответствии с определением надо доказать оценку(3.2.5). Заметим, чтоε0 = |u(x0 ) − y0 | = |u(x0 ) − u0 | = 0.47Глава 3. Задача Коши для обыкновенных дифференциальных уравненийИз оценки (3.2.6) имеемεi+1 6 (1 + c1 · h)εi + c2 · hk+1 66 (1 + c1 · h) · [(1 + c1 · h)εi−1 + c2 · hk+1 ] + c2 · hk+1 == (1 + c1 · h)2 · εi−1 + c2 · hk+1 · [1 + (1 + c1 · h)] 66 ...
66 (1 + c1 · h)i+1 · ε0 + c2 · hk+1 · [1 + (1 + c1 · h) + ... + (1 + c1 · h)i ] =i+1−1k+1 (1 + c1 · h)= c2 · h·.(1 + c1 · h) − 1Заметим теперь, что(1 + c1 · h)i+1 = exp{(i + 1) · ln(1 + c1 · h)} 6 exp{(i + 1) · c1 · h} =b−a· (i + 1)} 6 e c1 ·(b−a) .= exp{c1 ·nС помощью этого последнего неравенства из предыдущей цепочки неравенств следует оценкаεi+1[c2 · (e c1 ·(b−a) − 1)] k·h ,6c1выполняющаяся для всех i = 0, 1, ..., (n − 1). Сравнивая последнее неравенство с неравенством (3.2.5) и полагая C = [c2 · (e c1 ·(b−a) − 1)]/c1 , убеждаемсяв справедливости леммы.Теперь, в соответствии с леммой, чтобы доказать, что метод Эйлера является методом первого порядка точности, достаточно проверить неравенствоεi+1 6 (1 + c1 · h)εi + c2 · hk+1 ,при всех i = 0, 1, ..., n − 1.Вычитая (3.2.3) из равенств (3.2.4), получаемui+1 − yi+1h2 00= ui − yi + h · [f (xi , u(xi )) − f (xi , yi )] + u (ξi ).248(3.2.7)Глава 3.
Задача Коши для обыкновенных дифференциальных уравненийПо теореме Лагранжа о конечном приращении функции имеемf (xi , u(xi )) − f (xi , yi ) =∂f(xi , ηi ) · (ui − yi ),∂uгде точка ηi лежит между точками ui и yi .С помощью последнего равенства из соотношения (3.2.7) находим оценкуεi+1 6 εi + h · c1 · εi + c2 · h2 = εi · (1 + c1 h) + c2 · h2 ,где c1 = max |fu0 (x, u)|, c2 = 0.5 · max (|fx0 (x, u)| + |fu0 (x, u)|· | f (x, u)|)(x,u)∈D(x,u)∈D(см. формулу (3.2.2)).Таким образом, доказано, что явный метод Эйлера имеет первый порядокточности.3.2.2. Методы Рунге-КуттаРассмотрим теперь методы, погрешность которых при стремлении h кнулю убывает с более высокой скоростью.Метод Рунге-Кутта второго порядка точности. Его расчетные формулы:k1 = h · f (xi , yi ),k2 = h · f (xi + h, yi + k1 ),1yi+1 = yi + · (k1 + k2 ), i = 0, 1, ..., n − 1.249(3.2.8)Глава 3.
Задача Коши для обыкновенных дифференциальных уравненийМетод Рунге-Кутта четвертого порядка точности. Вычисления с помощью этого метода проводятся по формулам:k1 = h · f (xi , yi ),k1hk2 = h · f (xi + , yi + ),22hk2k3 = h · f (xi + , yi + ),22k4 = h · f (xi + h, yi + k3 ),yi+1 = yi +(3.2.9)1· (k1 + 2k2 + 2k3 + k4 ), i = 0, 1, ..., n − 1.6Мы рассмотрели методы, для которых при вычислении yi+1 нужно знатьлишь значение yi , а значения приближенного решения в предшествующихточках не входят в расчетные формулы. Иными словами, одношаговые методы — это методы с «короткой памятью».
Если «память метода» получше, тоего называют многошаговым. Более точно, метод численного интегрированиязадачи называется l-шаговым, если при вычислении значения yi+1 используются l величин yi−l+1 , yi−l+2 , ... , yi .3.2.3. Многошаговые методы АдамсаИз (3.1.1) следует, чтоZxi+1Zxi+1Zxi+1u(xi+1 ) − u(xi ) =u0 (x) dx =f (x, u(x)) dx ≈p(x) dx,xixixiгде p(x) — полином, аппроксимирующий f (x, u(x)). Пусть fi = f (xi , yi ), гдеyi — приближенное решение задачи (3.1.1), и в качестве p(x) возьмем интерполяционный полином, проходящий через l ранее найденных точек (xj , fj )(j = (i − l + 1), (i − l + 2), (i − l + 3), ... , i), включая текущую точку (xi , fi ).50Глава 3.
Задача Коши для обыкновенных дифференциальных уравненийЕсли l = 1, то имеем явный метод Эйлера (3.2.3). Если l = 2, то p(x) — линей-ûô ô é ôная функция (рис. 3.2,a), проходящая через две точки (xi−1 , fi−1 ) и (xi , fi ):yayáfi+1p(x)p(x)fifi-2fi+1fi-1fifi-1xi-20xi-1xixi+1x-hxi-10xihxi+12hxРис. 3.2.p(x) =(xi − x)(x − xi−1 )· fi−1 +· fi , p(xi−1 ) = fi−1 , p(xi ) = fi .hhИнтегрируя полином от xi до xi+1 , получаем двухшаговый метод Адамсавторого порядка точности (он также называется методом Адамса-Башфорта):yi+1 = yi +h· (3fi − fi−1 ).2(3.2.10)Если l = 3, то p(x) — парабола, проходящая через точки (xi−2 , fi−2 ),(xi−1 , fi−1 ) и (xi , fi ) (рис. 3.2,б ), а соответствующий трехшаговый методАдамса третьего порядка точности имеет видyi+1 = yi +h· (23fi − 16fi−1 + 5fi−2 ).12(3.2.11)Формулу (3.2.11) легко получить, если перейти к новой системе координат,51Глава 3.
Задача Коши для обыкновенных дифференциальных уравненийв которой координата x точки xi−1 равна нулю (см. рис. 3.2,б ). Тогдаp(x) = ax2 + bx + c;p(−h) = ah2 − bh + c = fi−2 ,p(0) = c = fi−1 ,p(h) = ah2 + bh + c = fi .Отсюда находимp(x) =fi−2 − 2fi−1 + fi 2 fi − fi−2·x +· x + fi−1 .2h22h(3.2.12)Интегрируя выражение (3.2.12) на отрезке [h, 2h], получим формулу(3.2.11).Если l = 4 то интерполяционный многочлен является кубическим и мыполучаем формулу Адамса четвертого порядка точностиyi+1 = yi +h· (55fi − 59fi−1 + 37fi−2 − 9fi−3 ).24(3.2.13)Многошаговые методы требуют в начале работы знания значений в первых l точках: y0 , y1 , ... , yl−1 . Мы не можем использовать, например, формулу(3.2.13) при i < 3.
Выход из положения состоит в применении какого-либо одношагового метода того же порядка точности, например метода Рунге-Кутта,до тех пор, пока не будет получено достаточное количество значений для проведения расчетов с помощью многошагового метода.Замечание. Наряду с рассмотренными явными методами существуют инеявные методы интегрирования дифференциальных уравнений. Приведемдва таких метода.Неявный метод Эйлера. Это метод первого порядка точности. Его расчетная формула:yi+1 = yi + h · f (xi+1 , yi+1 ), i = 0, 1, ..., n − 1.52(3.2.14)Глава 3. Задача Коши для обыкновенных дифференциальных уравненийЧтобы найти yi+1 надо решить это уравнение (может быть нелинейное) относительно этой переменной.Метод трапеций.
Это метод второго порядка точности. Значения yi+1находятся в результате решения уравненийyi+1 = yi +h· [f (xi , yi ) + f (xi+1 , yi+1 )], i = 0, 1, ..., n − 1.2(3.2.15)3.2.4. Правило Рунге практической оценки погрешностиЭто правило (см. гл. 1 и 2) применимо для практической оценки погрешности и при численном решении обыкновенных дифференциальных уравнений.(h)Пусть yi— значение в точке xi приближенного решения y (h) задачи Коши(3.1.1) на отрезке [a, b], найденное с шагом h, где xi = a + i · h, i = 0, 1, ..., n,(h/2)n — число разбиений отрезка [a, b], h = (b − a)/n.
И пусть yi— значениев той же точке xi , но приближенного решения y (h/2) , найденное с шагом h/2,т.е. число разбиений в этом случае равно 2n. Считается, что y (h/2) являетсярешением задачи Коши (3.1.1) с погрешностью ε, если(h/2)| yi(h)− yi |< ε,2k − 1где i = 1, 2, ..., n; k — порядок точности численного метода (например, k = 1— для метода Эйлера (3.2.3), k = 4 — для метода Рунге-Кутта (3.2.9)).Алгоритм вычислений. Допустим, что мы ищем численное решение задачи Коши (3.1.1) с помощью метода Рунге-Кутта четвертого порядка точности(k = 4). Опишем алгоритм вычислений, основанный на применении правилаРунге практической оценки погрешности.
Численное решение находят методом итераций. Пусть l — номер итерации, y l — численное решение, найденноес шагом hl , где hl — расчетный шаг на l-ой итерации. Очередную итерациюосуществляют следующим образом. Рассчитывают y l+1 с шагом hl+1 = hl /2.53Глава 3. Задача Коши для обыкновенных дифференциальных уравненийПосле этого проверяют выполнение неравенства| yil+1 − yil |< ε,2k − 1(3.2.16)строго говоря, во всех общих точках xi решений y l и y l+1 .Обычно выполнение неравенства (3.2.16) проверяют не во всех общих точках решений y l и y l+1 , а только в выделенных контрольных точках.
В качестве контрольных можно взять узловые точки {x0i }, соответствующие начальному числу разбиения n0 с шагом h0 : x0i = a + i · h0 , i = 0, 1, ..., n0 ,h0 = (b − a)/n0 . Число n0 (это целое число) определяется по формуле (см.(1.3.2) и (1.3.3))b−an0 = √+ 1.kε(3.2.17)Здесь k — порядок точности метода, квадратные скобки [ ], как и в (1.3.2) и(1.3.3), обозначают целую часть заключенного в них числа.h ib−a√√Замечание. Если k ε = b−ak ε , т.е. если дробная часть числаb−a√k εравняетсянулю, то 1 в формуле (3.2.17) можно не прибавлять.