Вычислительная математика лекции. Конспект лекций. Санкт петербург 2011 2 Оглавление
 Скачать 3.86 Mb.
|
|
10.2.3. Методы Рунге – Кутты. Приблизительно в 1900 г. немецкие математики Рунге и Кутта предложили способ построения формулы для φ(t,x,h), не содержащий производных от функции f(t,x). Суть предложенной модификации состоит в следующем. Требуется, используя формулу x k+1 = x k +hφ(t k ,x k ,h), разработать метод с порядком точности p. Будем искать φ(t k ,x k ,h) в виде 113 1 1 1 1 ( , , ) ; ( , ) ( , ); r=2, . p k k r r k k r r r k r k rs s s t x h c l l f t x l f t a h x h b l p Параметры c r , a r , b rs выбираются из условия совпадения максимально возможного числа слагаемых в разложении по степеням h для x(t k ) и x k +hφ(t k ,x k ,h). Использование методики поясним на примере p = 2. В этом случае φ(t,x,h) = c 1 f(t,x) +c 2 f(t + ha 2 , x+b 21 hf(t,x)). Разложим эту функцию в ряд по степеням h в окрестности точки (t k ,x(t k )). Учитывая, что x'(t k )=f(t k ,x k ), имеем k k 1 2- k k 2 2 2 21 t , x t , h c c f t , x t hc [a '( , ( )) '( , ( )) ( , ( ))] ( ) t k k x k k k k f t x t b hf t x t f t x t o h Ранее в процессе разложения решения x(t) в ряд Тэйлора при s = 2 было получено 2 k k k k 2 t , x t , h x'(t ) ''( ) ( ) 2 t , x t [ '( , ( )) 2 '( , ( )) ( , ( ))] ( ) k k t k k x k k k k h x t o h h f f t x t f t x t f t x t o h Чтобы обеспечить совпадение правых частей в двух последних выражениях, нужно выполнить следующие условия c 1 +c 2 =1, c 2 a 2 =1/2, c 2 b 21 =1/2.Считая c 2 произвольной константой, получаем семейство методов Рунге – Кутты второго порядка точности. Если c 2 =1/2, получаем метод Эйлера – Коши c 1 =c 2 =1/2, a 2 =b 21 =1 1 [ ( , ) ( , ( , ))]. 2 k k k k k k k k h x x f t x f t h x hf t x Другая распространенная разновидность метода получается при выборе c 2 =1 114 1 [ ( , ) ( , ( , ))]. 2 2 k k k k k k k k h h x x h f t x f t x f t x Широкое распространение получили методы Рунге – Кутты 4-ого порядка точности, а среди них следующая разновидность 1 1 2 3 4 1 2 1 3 2 4 3 ( 2 2 ) 6 ( , ) ( , ) 2 2 ( , ) 2 2 ( , ) k k k k k k k k k k h x x l l l l l f t x h h l f t x l h h l f t x l l f t h x hl На рисунке изображена область асимптотической устойчивости метода Рунге – Кутты четвертого порядка. 2 3 1 hIm λ hRe λ 115 10.2.4. Контроль точности методов Рунге – Кутты. В современных программах обязательно используются специальные алгоритмы автоматического выбора шага по условию обеспечения заданной точности. Один из распространенных подходов состоит в использовании метода Рунге ( метод двойного пересчета). Пусть в точке t k найдено решение x(t k ). Обозначим через 1 1 k h k x приближение к x(t k+1 ), найденное с помощью численного метода 1 ( , , ) (**) k k k k x x t x h h p – ого порядка точности. Ранее указывалось, что локальная погрешность на очередном шаге для методов Рунге – Кутты равна L k = α k 1 1 p k h . Следовательно, 1 1 1 1 1 ( ) k h p k k k k x t x h Уменьшим теперь шаг интегрирования вдвое и за два шага повторно вычислим приближенное решение в точке t k+1 . Обозначим это решение 1 /2 1 k h k x . Можно ожидать, что 1 1 /2 1 1 1 ( ) 2 2 k p h k k k k h x t x Используя этот и предыдущий результаты для исключения неизвестного α k , получаем 1 1 1 /2 /2 1 1 1 1 ( ) 2 1 k k k h h h k k k k p x x x t x Величину первого шага обычно задает пользователь. С помощью последней формулы оценивается величина погрешности. Если погрешность оказалась больше заданной, шаг уменьшают вдвое до тех пор, пока заданная погрешность будет достигнута. Если на очередном шаге погрешность оказалась существенно меньше заданной, шаг увеличивают в два раза. 116 10.3. Методы эквивалентных интегральных уравнений Метод эквивалентных интегральных уравнений основан на дуальности понятий дифференцирования и интегрирования: на идее эквивалентности обыкновенного дифференциального уравнения 0 0 ( ) ( , ), ( ) x t f t x x t x интегральному уравнению 0 t 0 t , ( ) x t x t f x d (1) Используя уравнение (1) (формулу Ньютона-Лейбница) на отрезке 1 [ , ] k k t t , получим формулу, которую можно считать базовой для построения последующих разностных схем 1 1 , ( ) k k t k k t x x f x d (2) При использовании для вычисления интеграла квадратурной формулы левых прямоугольников получим формулу явного метода ломаных Эйлера 1 ( , ) k k k k x x h f t x (3) При использовании квадратурной формулы правых прямоугольников получим формулу неявного метода ломаных Эйлера 1 1 1 ( , ) k k k k x x h f t x (4) При использовании квадратурной формулы трапеций получим формулу неявного метода трапеций 1 1 1 ( , ) ( , ) 2 k k k k k k h x x f t x f t x (5) Каковы дальнейшие пути повышения точности? На первый взгляд возможно использование более точных квадратурных формул, 117 например, формулы Симпсона, Гаусса и др. Однако все они требуют вычисления подынтегральной функции в некоторых точках внутри промежутка интегрирования, в то время как решение ( ) x t в этих точках еще не определено. Дальнейшее развитие метода эквивалентных интегральных уравнений связано с заменой подынтегральной функции интерполяционным полиномом, что приводит к методам Адамса. Джон Адамс (1819-1892г.) – английский математик и астроном. Основная идея методов Адамса состоит в замене подынтегральной функции на интервале 1 [ , ] k k t t интерполяционным полиномом. Первый подход состоит в использовании информации о ранее полученных точках 1 2 , , , k k k t t t Так используя результаты ранее выполненных шагов: узлы ( , ), k k t x 1 1 ( , ), k k t x …, ( , ) k m k m t x , а также значения подынтегральной функции в этих узлах: 1 , , , k k k m f f f , можно построить интерполяционный полином m-го порядка. При этом потребуется вычислить лишь одно значение функции ( , ) k k k f f t x . Для примера возьмем m=2 и, используя табличную функцию, куда внесем значения функции для 1 2 , , k k k t t t , 118 t t k t k-1 t k-2 f f k f k-1 f k-2 построим интерполяционный полином 2-го порядка 1 2 1 2 1 2 2 ( ) ( ) ( )( ) 2 k k k k k k k k k f f f f f p t f t t t t t t h h Подставив его под знак интеграла в (2), получим 1 2 1 2 1 2 0 1 2 ( ) 2 ( ) 2 23 16 5 12 12 12 k h k k t k k k t t t h k k k k k k k k k k x x p t dt x f f f f f f h d h h x h f f f Получили формулу метода 3-го порядка точности. Порядок точности определяется величиной s=m+1. Увеличивая число узлов и, соответственно, порядок m полинома, можно увеличивать порядок точности метода. Получающиеся этим методом соотношения принадлежат семейству многошаговых алгоритмов и относятся к явным методам Адамса, называемых также методами Адамса- Бэшфорта. Общий вид формул метода 1 0 m k k k r r x x h f mr s 1 2 3 4 m 0 1 2 3 β m0 1 3/2 23/12 55/24 β m1 -1/2 -16/12 -59/24 β m2 5/12 37/24 β m3 -9/24 119 В приведенной таблице указаны коэффициенты наиболее часто используемых формул. При использовании m=0 получаем формулу "явного метода ломаных Эйлера". Несомненным достоинством методов Адамса является тот факт, что все они независимо от своей точности требуют лишь однократного вычисления функции ( , ) f t x на одном шаге, и конкурировать с ними в этом плане трудно. Следует отметить, что область s=1 асимптотической устойчивости методов Адамса- Бэшфорта уменьшается с ростом s. (см. примерный рис.) Для вещественных отрицательных собственных значений ( ) i A ограничения на выбор шага имеют вид: ( ) 1.0 i h A − для метода Адамса второго порядка; ( ) 6 11 i h A − для метода третьего порядка; ( ) 0.3 i h A − для метода Адамса четвертого порядка. Первый подход при построении метода Адамса состоял в использовании информации о ранее полученных точках 1 2 , , , k k k t t t 2 Re , i h Im , i h 3 3 4 -- 2 120 Второй подход состоит в дополнении системы ранее вычисленных узлов ( , ), k k t x 1 1 ( , ), k k t x …, 1 1 ( , ) k m k m t x , а также значений подынтегральной функции в этих узлах: 1 1 , , , k k k m f f f дополнительным, подлежащим определению узлом 1 1 ( , ) k k t x и значением функции 1 k f для этого узла. Построив интерполяционный полином и подставив его в качестве подынтегральной функции в соотношение (2) получим формулы неявного метода Адамса или формулы Адамса-Мултона(Моултона) * 1 1 0 m k k k r r x x h f mr Коэффициенты формул Адамса-Мултона приведены в следующей таблице s 1 2 3 4 m 0 1 2 3 * 0 m 1 1/2 5/12 9/24 * 1 m 1/2 8/12 19/24 * 2 m -1/12 -5/24 * 3 m 1/24 Формула Адамса-Мултона первого порядка точности (при s=1) определяет уже известный нам неявный метод ломаных Эйлера, второго порядка точности (s=2) − совпадает с неявным методом трапеций (5). Методы Адамса-Мултона при 2 s являются А-устойчивыми, это значит, что условия асимптотической устойчивости численных решений линейных дифференциальных систем обеспечиваются с любым шагом при Re ( ) 0, i A i Многошаговые алгоритмы (в том числе явные методы Адамса при 2 s и неявные методы Адамса при 3 s ) не являются само стартующими, т.е. они требуют для начала интегрирования специальных стартовых 121 алгоритмов для расчета дополнительных начальных условий. В качестве стартовых методов может быть использован любой другой одношаговый метод, например один из методов Рунге- Кутты. Порядок сходимости стартового метода должен быть не ниже, чем у основного. 10.4 . Свойство жесткости систем дифференциальных уравнений Проблемы, связанные с решением "жестких" систем дифференциальных уравнений рассмотрим на примере решения двух асимптотически устойчивых задач Пример 1. Пусть матрица А размера 2 ×2 имеет 2 собственных значения 1 2 2, 1 Решение системы содержит две составляющие с экспонентами exp( 2 ) t и exp( ) t . Интервал решения определяется самой медленной экспонентой Т=1с. и обычно выбирается равным (3 5)T . Для построения графика решения вполне достаточно выбрать шаг временной дискретности 0.1 0.2 h и тогда график решения будет содержать десять-тридцать точек. Условие численной устойчивости при использовании явного метода ломаных Эйлера связано с выполнением неравенства max 2 i h , что полностью соответствует поставленной задаче. Пример 2. Пусть 4 1 2 10 , 1 . Решение вновь содержит две экспоненты и, т.к. самая медленная осталась прежней, время наблюдения решения остается сохраняется [0,3] t . Но на этот раз на начальном участке имеется резко меняющаяся экспонента с малой постоянной времени 4 1 1 1 10 с . Для построения графика такого решения необходимо взять десяток точек с шагом 5 10 h , а затем, на оставшейся, большей части графика, для решения использовать шаг порядка нескольких десятых. Однако условие устойчивости требует выбора шага |