авторефераты диссертаций БЕСПЛАТНАЯ БИБЛИОТЕКА РОССИИ

КОНФЕРЕНЦИИ, КНИГИ, ПОСОБИЯ, НАУЧНЫЕ ИЗДАНИЯ

<< ГЛАВНАЯ
АГРОИНЖЕНЕРИЯ
АСТРОНОМИЯ
БЕЗОПАСНОСТЬ
БИОЛОГИЯ
ЗЕМЛЯ
ИНФОРМАТИКА
ИСКУССТВОВЕДЕНИЕ
ИСТОРИЯ
КУЛЬТУРОЛОГИЯ
МАШИНОСТРОЕНИЕ
МЕДИЦИНА
МЕТАЛЛУРГИЯ
МЕХАНИКА
ПЕДАГОГИКА
ПОЛИТИКА
ПРИБОРОСТРОЕНИЕ
ПРОДОВОЛЬСТВИЕ
ПСИХОЛОГИЯ
РАДИОТЕХНИКА
СЕЛЬСКОЕ ХОЗЯЙСТВО
СОЦИОЛОГИЯ
СТРОИТЕЛЬСТВО
ТЕХНИЧЕСКИЕ НАУКИ
ТРАНСПОРТ
ФАРМАЦЕВТИКА
ФИЗИКА
ФИЗИОЛОГИЯ
ФИЛОЛОГИЯ
ФИЛОСОФИЯ
ХИМИЯ
ЭКОНОМИКА
ЭЛЕКТРОТЕХНИКА
ЭНЕРГЕТИКА
ЮРИСПРУДЕНЦИЯ
ЯЗЫКОЗНАНИЕ
РАЗНОЕ
КОНТАКТЫ


Pages:     | 1 || 3 | 4 |   ...   | 9 |

«Е.А. Новиков, Ю.В. Шорников КОМПЬЮТЕРНОЕ МОДЕЛИРОВАНИЕ ЖЕСТКИХ ГИБРИДНЫХ СИСТЕМ Е.А. Новиков, Ю.В. Шорников КОМПЬЮТЕРНОЕ ...»

-- [ Страница 2 ] --

2.2. КОНТРОЛЬ ТОЧНОСТИ ВЫЧИСЛЕНИЙ Здесь и ниже будут рассматриваться оценки локальной и глобаль ной ошибок в смысле главного члена, т. е. далее будет оцениваться только первый член при разложении соответствующих ошибок в ряды Тейлора по степеням h. Итак, пусть для решения задачи (2.2) имеется метод вида (2.5) p -го порядка точности и пусть его локальная ошибка n, p представима в виде n, p = h p +1(tn ) + O(h p+ 2 ), (2.10) где – некоторая функция, не зависящая от размера шага интегриро вания.

Для записи (2.10) требуется определенная гладкость правой части исходной задачи на промежутке [tn, tn +1 ], которая предполагается вы полненной. Известно, что при условии (2.10) глобальную ошибку n, p можно вычислить по формуле Г л а в а 2. КОНТРОЛЬ ТОЧНОСТИ И УСТОЙЧИВОСТИ ОДНОШАГОВЫХ МЕТОДОВ n, p = h p x(tn ) + O(h p+1 ), (2.11) где x(t ) есть решение задачи Коши x(t ) f (t, y (t )) x(t ) = (t ), x(t0 ) = 0, (2.12) f = f / y – матрица Якоби системы дифференциальных уравнений.

Непосредственное использование оценки (2.11) в неравенстве для контроля точности вычислений может приводить к значительному увеличению вычислительных затрат. Дело в том, что для определения n, p требуется оценка локальной ошибки (2.10), что означает, как пра вило, дополнительные вычисления правой части исходной задачи.

Кроме того, решение задачи (2.12) связано с вычислением N 2 компо нент матрицы Якоби системы (2.2) плюс затраты на интегрирование (2.12). Поэтому такой способ используется лишь иногда при решении жестких задач, когда матрица Якоби уже известна.

Обычно при реализации A -устойчивых или L -устойчивых мето дов вида (2.5) имеется матрица Dn вида Dn = E ahf (tn, yn ), (2.13) где E – единичная матрица, a – известное число.

Поэтому ниже будем считать матрицу Dn заданной. Опишем один способ оценки глобальной ошибки методов с неограниченной об ластью устойчивости, который в ряде случаев приводит к незначи тельному увеличению вычислительных затрат. Для этого докажем сле дующую теорему.

Теорема 2.1. Пусть для решения задачи (2.1) имеется метод p -го порядка точности и пусть его локальная ошибка n, p представима в виде n, p = h p +1 f (tn, y (tn )) (tn ) + O(h p + 2 ), (2.14) где (t ) есть некоторая функция, не зависящая от размера шага ин тегрирования. Тогда для любого номера шага n имеет место соот ношение || x(tn ) xn || = O(h), (2.15) 2.2. Контроль точности вычислений где || || – некоторая норма в R N, x(t ) – решение задачи Коши (2.12), а xn вычисляется по следующей рекуррентной формуле:

xn +1 = xn + a 1 ( n xn ) + Dn 1 ( xn n ), (2.16) x0 = 0, n = (tn ).

Для доказательства подставим (t ) = f (t, y (t )) (t ) в (2.12), т. е. для определения x(t ) получим следующую задачу:

x(t ) = f (t, y (t ))[ x(t ) (t )], x(t0 ) = 0. (2.17) Численное решение (2.17) определим с помощью схемы типа Розен брока первого порядка точности, которая применительно к задаче (2.1) имеет вид Dn ( yn+1 yn ) = hf (tn, yn ), (2.18) где величина a при задании Dn по формуле (2.13) определена при описании исходного метода p -го порядка точности. Применяя (2.18) для решения (2.17), запишем Dn ( xn +1 xn ) = hf (tn, yn ) xn hf (tn, yn ) n. (2.19) Добавляя в правую часть равенства (2.19) соотношение a ( n n + xn xn ) и группируя соответствующим образом слагае мые, получим Dn ( xn+1 xn ) = a 1[( xn n ) + Dn ( n xn )]. Применяя слева к полученному соотношению матрицу Dn 1, получим формулу (2.16). Равенство (2.15) следует из сходимости схемы (2.18) с первым порядком. Теорема доказана.

Заметим, что схема (2.18) при a 0,5 является A -устойчивой, а при a = 1 L -устойчивой. Это означает, что при значениях a, 0 a 0,5, оценка глобальной ошибки, вычисленная с применением (2.16), может возрастать, в то время как истинная ошибка будет убы вать. Однако для многих L -устойчивых методов (2.5) имеет место 0 a 0,5, и поэтому для них использование (2.16) при выборе вели чины шага интегрирования приведет к понижению эффективности Г л а в а 2. КОНТРОЛЬ ТОЧНОСТИ И УСТОЙЧИВОСТИ ОДНОШАГОВЫХ МЕТОДОВ алгоритма интегрирования. С целью устранения этого недостатка рас смотрим (2, 1)-метод, который применительно к (2.1) имеет вид yn+1 = yn + ak1 + (1 a )k2, Dn k1 = hf (tn, yn ), Dn k2 = k1. (2.20) Нетрудно убедиться, что при 0,5 0, 25 2 a 0,5 + 0, 25 2 данная схема имеет первый порядок точности и является L -устойчивой.

Теорема 2.2. Пусть выполнены условия теоремы 2.1. Тогда имеет место соотношение (2.15), где 2a 1 1 a 1 xn +1 = n + Dn 1 ( xn n ) + Dn ( xn n ), a a (2.21) x0 = 0, n = (tn ).

Для доказательства применим численную схему (2.20) для решения задачи (2.17). Получим Dn ( xn +1 xn ) = ahf (tn, yn )[ xn n ] + + (1 a) Dn 1hf (tn, yn )[ xn n ].

Добавим в правую часть соотношение ( n xn n + xn ). Учиты вая вид матрицы Dn, сгруппируем определенным образом слагаемые и умножим слева полученное соотношение на Dn. В результате запишем Dn ( xn +1 n ) = Dn ( xn n ) + a 1 (1 a ) [ ahf (tn, yn ) xn ahf (tn, yn ) n ].

Теперь добавим к выражению в квадратных скобках соотношение ( n xn n + xn ) и снова сгруппируем слагаемые, учитывая вид Dn.

Получим Dn ( xn +1 n ) = a 1 (2a 1) Dn ( xn n ) + a 1 (1 a)( xn n ).

Умножив слева полученное уравнение на матрицу Dn 2, имеем соотношение (2.21). Равенство (2.15) следует из сходимости (2.20).

Теорема доказана.

Теперь рассмотрим вопрос о том, в каком случае применение (2.16) для оценки глобальной ошибки более выгодно с точки зрения вычис 2.2. Контроль точности вычислений лительных затрат, чем непосредственное интегрирование (2.12). Преж де всего отметим, что при вычислении стадий A -устойчивых или L устойчивых методов (2.5) будет использоваться LU -разложение мат рицы Dn с выбором главного элемента по строке или столбцу, а ино гда по всей матрице, т. е. при вычислении приближенного решения по формуле (2.5) осуществляется декомпозиция матрицы Dn. Поэтому вычисление оценки ошибки по формуле (2.16) не приводит к значи тельному увеличению вычислительных затрат, которые практически определяются числом операций на выполнение обратного хода в мето де Гаусса. Далее, как уже отмечалось выше, использование (2.12) свя зано с необходимостью оценки локальной ошибки, что, как правило, приводит к дополнительным вычислениям правой части задачи (2.2) и как следствие – увеличению вычислительных затрат. При использова нии (2.16) требуется оценить функцию (t ), что в ряде случаев удает ся осуществить бесплатно. Ниже при оценке (t ) будем поступать следующим образом.

Пусть для решения задачи (2.2) имеются два метода вида (2.5) ( p 1) -го и p -го порядков точности. Обозначим через yn, p 1 и yn, p соответствующие приближения к решению, а через n, p 1 и n, p – их локальные ошибки. Выберем коэффициенты методов такими, чтобы локальные ошибки были согласованы:

n, p 1 = c p 1h p p (tn ) + O (h p +1 ), n, p = = c p h p +1 f (tn, y (tn )) p (tn ) + O (h p + 2 ), (2.22) где c p 1 и c p – некоторые вычисляемые постоянные, а функция p (t ) не зависит от размера шага интегрирования. Вычисления пред полагается проводить по методу p -го порядка, для определения гло бальной ошибки которого требуется оценка функции (t ) = c p p (t ).

Отметим, что требование (2.22) принципиально возможно, потому что в разложении ошибки в ряд Тейлора можно оставить только те члены, которые получаются на линейной задаче. Ясно также, что (2.22) означает некоторые дополнительные соотношения на коэффициенты схемы (2.5). Однако для выполнения условий (2.22) во многих случаях Г л а в а 2. КОНТРОЛЬ ТОЧНОСТИ И УСТОЙЧИВОСТИ ОДНОШАГОВЫХ МЕТОДОВ могут быть использованы те коэффициенты, которые остаются сво бодными после удовлетворения требованиям аппроксимации. Класс методов (2.5) с условиями (2.22) не пуст. Теперь докажем утвержде ние, с использованием которого далее будет определяться оценка функции p (t ).

Теорема 2.3. Пусть для решения задачи (2.2) имеется два метода вида (2.5) ( p 1) -го и p -го порядков точности и пусть их локальные ошибки удовлетворяют требованиям (2.22). Тогда выполняется соот ношение p (tn ) c 11h p ( yn, p yn, p 1 ) = O(h). (2.23) p Для доказательства предположим, что приближение к решению yn в точке tn вычислено методом p -го порядка, т. е. выполнено равенство yn = y (tn ) + n, n = O(h p ). (2.24) Тогда приближенное решение yn+1 в точке tn +1, вычисленное методом p -го порядка, можно представить в виде yn+1, p = yn + h f (tn, yn, h) = y (tn ) + n + h f (tn, y (tn ) + n, h).

Разлагая функцию f (tn, y (tn ) + n, h) в ряд Тейлора в окрестности точки y (tn ) и используя очевидное равенство f (tn, y (tn ), h) / y f (tn, y (tn )) / y = O(h), (2.25) получим yn+1, p = y (tn ) + n + h f (tn, y (tn ), h) + hf (tn, y (tn )) n + O (h p + 2 ).

Добавляя разность y (tn+1 ) y (tn +1 ) в правую часть полученного соот ношения и учитывая определение погрешности аппроксимации, будем иметь yn+1, p = y (tn +1 ) n +1, p + n + hf (tn, y (tn )) n + O(h p + 2 ). (2.26) 2.2. Контроль точности вычислений Проводя аналогичные рассуждения для метода ( p 1) -го порядка, имеем yn+1, p 1 = y (tn+1 ) n+1, p 1 + n + hf (tn, y (tn )) n + O(h p + 2 ). (2.27) Учитывая, что n, p = O(h p +1 ), n, p 1 = O (h p ), и вычитая (2.27) из (2.26), получим n +1, p 1 = yn +1, p yn +1, p 1 + O (h p +1 ), (2.28) откуда, с учетом (2.22), следует (2.23). Теорема доказана.

При использовании описанного выше способа оценки глобальной ошибки экономия вычислительных затрат достигается за счет того, что в ряде случаев при построении метода ( p 1) -го порядка точности ис пользуются те же самые стадии, что и в методе p -го порядка, т. е. ме тод ( p 1) -го порядка является побочным продуктом метода p -го по рядка.

Оценка локальной ошибки с помощью вложенных методов осуще ствлялась рядом авторов. При этом численная формула более высокого порядка точности была вспомогательным методом и он строился для оценки локальной ошибки метода более низкого прядка. Поэтому та кой подход всегда приводил к увеличению вычислительных затрат.

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

Возможны и другие способы оценки функции p (t ), не связанные с построением вспомогательного метода. При построении конкретных методов этот вопрос будет обсуждаться более подробно. Вычисление оценки глобальной ошибки n, p метода вида (2.5) p -го порядка точ ности по формуле n, p = h p xn, (2.29) где xn определяется из соотношения (2.16), связано с вычислением матрицы Якоби. Во многих случаях при решении задач умеренной же сткости можно применить явные методы, в которых матрица Якоби не Г л а в а 2. КОНТРОЛЬ ТОЧНОСТИ И УСТОЙЧИВОСТИ ОДНОШАГОВЫХ МЕТОДОВ используется и для которых оценка ошибки по формуле (2.16) означает существенное увеличение вычислительных затрат. В то же время при интегрировании умеренно жестких задач явными методами основные затраты приходятся на участок установления, где шаг интегрирования ограничен не требованием точности, а устойчивостью численной схе мы и где точность расчетов, как правило, значительно завышена. Как показывают расчеты, в этом случае в качестве оценки ошибки метода можно использовать величину c p h p p (t ). Докажем следующее ут верждение.

Теорема 2.4. Пусть для решения скалярной задачи (2.1) имеется метод p -го порядка точности, локальная ошибка которого имеет вид (2.22). Тогда выполняется неравенство || x(t ) || | c p | || p (t ) ||, где || x(t ) || = max t0 z t | x( z ) | и || p (t ) || = max t0 z t | p ( z ) |, а x(t ) есть решение задачи Коши x(t ) f (t, y (t )) x(t ) = c p f (t, y (t )) p (t ), x(t0 ) = 0.

t f ( z, y ( z ))dz. В силу Для доказательства введем обозначение F = t знакопостоянства f (t, y (t )) для скалярного уравнения вида (2.1) с ис пользованием теоремы о среднем можно записать t || x(t ) || = || c p e F e F f ( z, y ( z ) ) p ( z )dz || t t | c p | | p ( z0 ) | || e F e F dF || | c p | || p (t ) ||, t где z0 [t0, tk ]. Теорема доказана.

Основываясь на этом соображении при N 1, будем предполагать выполнение неравенства max1i N || xi (t ) || | c p | max1i N || pi (t ) ||, т. е. в качестве оценки глобальной ошибки n, p используем иногда 2.2. Контроль точности вычислений следующую оценку: n, p c p h p p (tn ). Учитывая теорему 2.3, эту оценку можно записать в виде n, p c p c 11 ( yn, p yn, p 1 ). (2.30) p Теперь для контроля точности вычислений и при выборе величины шага интегрирования можно применить неравенство || n, p ||, (2.31) где – требуемая точность интегрирования, || || – некоторая норма в R N. Во всех приведенных ниже алгоритмах интегрирования левая часть неравенства (2.31) определяется по формуле || n, p || = max1i N | i, p | /(| yn | + r ), где i – номер компоненты, r – i n положительный параметр. Если по i -й компоненте решения выполня i ется неравенство | yn | r, то контролируется абсолютная ошибка r, в противном случае – относительная ошибка.

Отметим еще одну важную особенность построенных оценок. Для этого рассмотрим L -устойчивый метод (2.5), который применим для решения задачи (2.7). Тогда в (2.8) в силу L -устойчивости схемы вы полняется Q ( z ) 0 при z. Так как для точного решения y (tn+1 ) = exp(h) y (tn ) задачи (2.7) выполняется аналогичное свойство, естественно требование стремления к нулю глобальной ошибки при условии z. В силу того что все приведенные выше рассуждения относились к оценке главного члена ошибки, это требование может не выполняться. Поэтому с целью исправления асимптотического пове дения глобальной ошибки вместо построенных выше оценок n, p будем рассматривать оценку n, p ( jn ) вида n, p ( jn ) = D1 jn n, p, n 1 jn J f, где параметр J f зависит от асимптотических свойств n, p. Теперь для контроля точности и при выборе величины шага ин тегрирования можно проверять следующее неравенство:

|| D1 jn n, p ||, 1 jn J f. (2.32) n Г л а в а 2. КОНТРОЛЬ ТОЧНОСТИ И УСТОЙЧИВОСТИ ОДНОШАГОВЫХ МЕТОДОВ Используя (2.13), нетрудно видеть, что в смысле главного члена, т. е. в смысле первого члена при разложении ошибок в ряды Тейлора по степеням h, n, p и n, p ( jn ) совпадают при любом значении jn, 1 jn J f. Для явных методов, в которых матрица Якоби не участву ет в вычислительном процессе, jn 1. Отметим, что применение (2.32) вместо (2.31) к существенному увеличению вычислительных затрат не приводит. При условии z 0 оценка n, p (1) = n, p правильно отража ет поведение глобальной ошибки и нет смысла проверять (2.32) при других значениях jn. При резком увеличении шага поведение ошибки n, p может оказаться неудовлетворительным, что проявляется в неоп равданном уменьшении величины шага и повторных вычислениях ре шения (возвратах). Проверка неравенства (2.32) при jn 1 осуществ ляется достаточно редко. Как показывают расчеты, использование (2.32) вместо (2.31) позволяет сократить вычислительные затраты при мерно на 10...20 %. При практической реализации алгоритмов интег рирования неравенство (2.32) используется так. При каждом фиксиро ванном n выбирается наименьшее значение jn, при котором выполняется неравенство (2.32). Если оно не выполняется ни при ка ком значении jn, включая jn = J f, то шаг уменьшается и происходит повторное вычисление решения.

Теперь приведем формулу для выбора величины шага интегриро вания с использованием неравенства (2.32). Пусть приближение к ре шению yn вычислено с шагом hn. Учитывая, что n, p ( jn ) = O(h p ), определим параметр hnew по формуле hnew = qhn, где значение q вы числяется из уравнения q p || n, p ( jn ) || =. Если q 1, т. е. неравенство (2.32) нарушается при всех jn, 1 jn J f, и, следовательно, требуе мая точность не выполняется, то hn полагается равным hnew и проис ходит повторное вычисление решения. Если q 1, т. е. при некотором значении jn неравенство (2.32) выполнено и, следовательно, требуе мая точность достигнута, то hn +1 полагается равным hnew и выполня ется следующий шаг интегрирования.

2.3. Контроль устойчивости 2.3. КОНТРОЛЬ УСТОЙЧИВОСТИ При численном решении задачи Коши для жестких систем обыкно венных дифференциальных уравнений, как правило, используются не явные или полуявные численные схемы. При реализации этих формул на ЭВМ требуется вычисление и обращение матрицы Якоби. В случае большой размерности исходной задачи (2.1) возникают трудности как с размещением элементов этой матрицы в оперативной памяти ЭВМ, так и с ее обращением.

Этого недостатка лишены явные формулы, но они в большинстве своем не пригодны для решения жестких задач по следующим причи нам. Обычно алгоритм управления величиной шага строится на кон троле точности численной схемы. Однако при решении жестких задач явными методами на большой части интервала интегрирования шаг ограничен не точностью, а устойчивостью численной формулы. В ре зультате на этой части интервала происходит раскачивание величины шага. Из-за отсутствия ограничения на шаг сверху он может увели читься настолько сильно, что рост ошибок из-за возникающей неус тойчивости приводит к аварийной остановке вследствие переполнения арифметического устройства. В лучшем случае почти каждый шаг со провождается повторным вычислением решения, что приводит к поте ре эффективности алгоритма интегрирования. Этого можно избежать, если наряду с точностью контролировать и устойчивость численной схемы. Отметим, что здесь не ставится вопрос о замене A -устойчивых или L -устойчивых методов явными. Речь идет лишь о тех умеренно жестких задачах, где явные схемы эффективнее методов с неограни ченной областью устойчивости.

Из результатов расчетов жестких дифференциальных уравнений алгоритмами на основе явных формул с выбором шага по точности и устойчивости следует, что на тех участках интервала интегрирования, где устойчивость оказывает решающую роль на размер шага, точность вычислений получается значительно выше требуемой, хотя вычисли тельные затраты меньше, чем при интегрировании алгоритмом без контроля устойчивости. Это естественно, потому что локальная ошиб ка убывает вследствие убывания производных решения, а предыдущая ошибка уменьшается за счет контроля устойчивости. В такой ситуации предпочтительнее считать по численной формуле более низкого по рядка точности, если ее область устойчивости существенно шире Г л а в а 2. КОНТРОЛЬ ТОЧНОСТИ И УСТОЙЧИВОСТИ ОДНОШАГОВЫХ МЕТОДОВ области устойчивости метода более высокого порядка. В качестве кри терия переключения с одной формулы на другую можно использовать неравенство для контроля устойчивости.

Еще одним соображением в пользу контроля устойчивости являет ся следующее. При интегрировании задачи Коши для жестких систем дифференциальных уравнений алгоритмами на основе неявных или полуявных формул основную часть в общих вычислительных затратах занимает время декомпозиции матрицы Якоби. Жесткость задачи мо жет быть настолько велика, что использование явной формулы даже с контролем устойчивости приводит к столь гигантскому числу шагов интегрирования, что такая задача становится непосильной для совре менных ЭВМ. Однако там, где устойчивость не влияет на размер шага, более естественно применять алгоритмы на основе явных формул. Там же, где условие устойчивости накладывает жесткие ограничения на шаг, целесообразно использовать полуявные или неявные схемы. Та ким образом, эффективность алгоритма интегрирования можно повы сить за счет расчетов переходных участков по явной формуле, которая не предполагает вычисление и обращение матрицы Якоби. В качестве критерия перехода с одной формулы на другую естественно использо вать неравенство для контроля устойчивости.

Теперь перейдем к получению неравенства для контроля устойчи вости методов (2.5). Применяя (2.5) к (2.7), получим yn+1 = Q( z ) yn, z = h. Условие устойчивости (2.5) имеет вид | Q( z ) | 1, кривая | Q( z ) | = 1 ограничивает область устойчивости чис ленной схемы в плоскости z = h. Учитывая, что при всех h из дан ной области метод устойчив, условие устойчивости можно еще запи сать в виде | h | D, где положительная постоянная D связана с размером области устойчивости. В случае, когда размерность задачи N 1, условие устойчивости имеет вид | i h | D, где i, 1 i N, есть собственные числа матрицы Якоби.

В настоящее время можно выделить два подхода к контролю ус тойчивости. Первый способ связан с оценкой максимального собст венного числа матрицы Якоби f (t, y (t ), через ее норму с последую щим контролем неравенства h || f (t, y (t ) || D. Ясно, что для явных методов, где матрица Якоби не участвует в вычислительном процессе, это приводит дополнительно к ее нахождению и, следовательно, к уве личению вычислительных затрат.

2.3. Контроль устойчивости Второй подход основан на оценке максимального собственного числа max матрицы Якоби степенным методом через приращения правой части исходной задачи с последующим контролем неравенства h | max | D. Во всех рассмотренных ситуациях этот подход не при водит к значительному увеличению вычислительных затрат. Изложим этот способ подробно. Вычитая (2.5) из (2.6), получим n+1 = n + h[ f (tn, y (tn ), h) f (tn, yn, h)] + n +1, (2.33) где n = y (tn ) yn. Пусть размерность N задачи равна единице и пусть к функции f применима теорема о среднем относительно y.

Обозначим Q(tn, xn, h) = E + h f (tn, xn, h) / y. Тогда из формулы (2.33) имеем n+1 = Q (tn, xn, h) n + n +1, xn [ y (tn ), yn ]. (2.34) В случае произвольного N вместо (2.34) можно записать || n +1 || || Q(tn, xn, h) n || + || n +1 ||, (2.35) где || || – некоторая норма в R N. Из (2.35) следует, что величина гло бальной ошибки n+1 в точке tn+1 зависит от величины локальной ошибки n+1 на данном шаге и от величины преобразованной преды дущей ошибки n. Для того чтобы схема (2.5) была устойчивой, дос таточно, чтобы оператор Q (tn, xn, h) не усиливал предыдущую ошиб ку. Ясно, что это требование будет выполнено, если при каждом n контролировать неравенство || Q(tn, xn, h) n || D1 || n ||, (2.36) где D1 – постоянная, 0 D1 1. Непосредственное использование (2.36) для контроля устойчивости связано с вычислением оператора Q (tn, xn, h) в некоторой неизвестной точке xn, что может приводить к значительным вычислительным затратам. Поступим следующим обра зом.

Пусть снова для решения задачи (2.1) имеется два метода вида (2.5) ( p 1) -го и p -го порядков точности и пусть их локальные ошибки согласованы:

Г л а в а 2. КОНТРОЛЬ ТОЧНОСТИ И УСТОЙЧИВОСТИ ОДНОШАГОВЫХ МЕТОДОВ n, p 1 = c p 1h p f ( p 1) f + O( h p + 2 ), (2.37) n, p = c p h p +1 f p f + O(h p + 2 ), где элементарные дифференциалы вычислены в точке (tn, y (tn )). Заме тим, что условие (2.37) более жесткое, чем (2.22), и поэтому для мето дов вида (2.5) с требованием (2.37) все рассуждения относительно кон троля точности вычислений остаются справедливыми, причем p = f ( p 1) f.

Величину максимального собственного числа c max при условии n (2.37) можно оценить с помощью степенного метода, т. е.

{ } {c p in, p1}, 1 i N.

h c max c p 1i, p (2.38) n n Тогда для контроля устойчивости численной формулы (2.5) можно ис пользовать, например, следующее неравенство:

max | i, p | / | in, p 1 | | c p / c p 1 | D, (2.39) n 1i N где постоянная D связана с размером области устойчивости. Однако для получения (2.39) требуется выполнение (2.37), что приводит к до полнительным соотношениям на коэффициенты численной схемы и, возможно, к дополнительным вычислениям правой части исходной задачи. Дело в том, что если для контроля точности вычислений требо валась лишь оценка локальной ошибки метода ( p 1) -го порядка, то для получения неравенства для контроля устойчивости дополнительно требуется оценка локальной ошибки метода p -го порядка. Это, как правило, приводит к увеличению вычислительных затрат.

С другой стороны, при получении неравенства для контроля устой чивости параметр p в условиях (2.37) определяет только число итера ций в степенном методе. Поэтому требования типа (2.37) можно рас сматривать при некотором значении l, l p. За счет этого можно избежать некоторых дополнительных вычислений правой части ис ходной задачи (2.1). При построении конкретных алгоритмов интегри рования этот вопрос будет рассмотрен более подробно. Там неравенст во типа (2.39) будет получено несколько иначе, хотя идейная сторона 2.3. Контроль устойчивости вопроса останется прежней оценка максимального собственного числа степенным методом.

Отметим следующее важное обстоятельство. Для контроля устой чивости можно использовать и другие неравенства, отличающиеся от (2.39). Например, можно рассмотреть следующее:

| n, p | / max | i, p 1 | | c p / c p 1 | D, j (2.40) n 1i N где j есть значение i, при котором достигается максимум знаменате ля. Использование того или иного неравенства зависит от цели, кото рой нужно достичь. Если нужно, чтобы устойчивость играла бльшую роль, то используется неравенство (2.39). Если же нужно построить алгоритм, который работал бы не хуже, чем без контроля устойчиво сти, то целесообразнее использовать неравенство (2.40), так как в этом случае, по крайней мере для симметричных матриц Якоби, модуль максимального собственного числа, вычисленный с помощью (2.40), будет меньше истинного.

Оценка (2.38) является грубой. Это связано с тем, что максималь ное собственное число матрицы Якоби вовсе не обязательно сильно отделено от остальных. В то же время нельзя без увеличения вычисли тельных затрат значительно увеличить число итераций в степенном методе, так как это связано с необходимостью получения условий типа (2.37) при достаточно большом значении параметра p. Кроме того, как видно из (2.37), дополнительные искажения в оценку максималь ного собственного числа будет вносить нелинейность задачи (2.1). По n n max = i =1 n max n, где этому в ряде случаев более надежна оценка cp c величина c max, 1 i n, определяется по формуле (2.38);

n – номер n текущей точки интегрирования. Кроме того, в алгоритме интегрирова ния можно предусмотреть обращение к подпрограмме, в которой оценку задает вычислитель. Ниже оценку максимального собственного числа матрицы Якоби, вычисленную одним из приведенных выше спо собов, будем обозначать n max. Для контроля устойчивости (2.5) можно использовать неравенство | hn n max | D, (2.41) где hn есть текущий шаг интегрирования.

Г л а в а 2. КОНТРОЛЬ ТОЧНОСТИ И УСТОЙЧИВОСТИ ОДНОШАГОВЫХ МЕТОДОВ 2.4. РЕАЛИЗАЦИЯ ЯВНЫХ МЕТОДОВ Далее речь пойдет о применении явных методов для решения задач умеренной жесткости. В явных методах не используется матрица Яко би, и поэтому прибегнем к оценке ошибки типа (2.30). Отметим, что при получении (2.30) может не применяться информация в точке tn +1, потому что вовсе не обязательно, чтобы даже одна стадия метода вы числялась в точке (tn +1, yn +1 ). При решении жестких задач в случае быстрого изменения решения это может стать причиной неудовлетво рительной точности. Поэтому при получении оценки (2.30) следует привлекать f (tn +1, yn +1 ). В случае успешного шага интегрирования это не приведет к существенному увеличению вычислительных затрат, потому что f (tn +1, yn +1 ) используется при выполнении следующего шага интегрирования.

Ниже будем предполагать, что имеются две оценки глобальной ошибки:, p и, p, причем при вычислении, p используется n n n f (tn +1, yn +1 ). В принципе, по вычислительным затратам оценка, p n, p, но для некоторых методов они могут совпадать.

дешевле, чем n Теперь для контроля точности вычислений и при выборе величины шага будем контролировать неравенства || n, p ||, || n, p ||, (2.42) где || || – некоторая норма в R N, – требуемая точность расчетов.

Первое неравенство (2.42) можно использовать для контроля точности вычислений, а при выборе величины шага дополнительно проверять второе неравенство. Более жесткий контроль при выборе шага, чем при контроле точности, позволит избежать некоторых повторных вычисле ний решения (возвратов) вследствие невыполнения точности. Иногда для повышения надежности контроля точности имеет смысл проверять оба неравенства (2.42).

Обозначим через vn оценку величины h | n max |, где n max есть максимальное собственное число матрицы Якоби системы (2.1). Тогда 2.4. Реализация явных методов vn D, где D связано для контроля устойчивости имеем неравенство с размером области устойчивости метода.

Пусть приближение к решению yn в точке tn вычислено с шагом hn. Это означает, что выполнено первое неравенство (2.42). Если оп ределить параметр q1 по формуле q1p || n, p || =, то q1 1, где p – порядок точности метода. Теперь новый шаг hnew1 из условия точно сти можно вычислить по формуле hnew1 = min(q1, q2 )hn, где q2 опре p деляется из уравнения q2 || n, p || =. Даже в случае успешного вы числения приближения к решению новый шаг интегрирования может быть уменьшен, потому что величина q2 может быть меньше едини цы. В этом случае не имеет смысла контролировать устойчивость, а приближение yn+1 в точке tn +1 вычисляется с шагом hn +1 = hnew1.

Величину шага hnew2 из условия устойчивости определим по фор муле hnew2 = rhn, где параметр r с учетом равенства vn = O(h) вы числяется из уравнения rvn = D. Предполагая, что поведение решения не сильно меняется от шага к шагу, величину прогнозируемого шага hn +1 можно вычислить по формуле hn +1 = min(hnew1, hnew2 ) / 1,1. По стоянная 1,1 внесена для того, чтобы создать некоторый запас по точ ности и устойчивости и тем самым избежать некоторых повторных вычислений решения.

Итак, если неравенства (2.42) надежны, а спектр матрицы Якоби и производные решения от шага к шагу изменяются незначительно, то приближение к решению yn +1 в точке tn +1 с шагом hn +1 должно быть вычислено успешно. Так как оценка максимального собственного чис ла является грубой, то это может приводить к неоправданному колеба нию величины шага, что не может не сказываться на эффективности алгоритма интегрирования. Поэтому в практических алгоритмах ин тегрирования контроль устойчивости используется как некоторый ог раничитель на рост шага, т. е. при q2 1 величина прогнозируемого шага вычисляется по формуле hn +1 = max( hn, min[q1, q2, r ]hn ). (2.43) Г л а в а 2. КОНТРОЛЬ ТОЧНОСТИ И УСТОЙЧИВОСТИ ОДНОШАГОВЫХ МЕТОДОВ Эта формула более инерционна, что позволяет несколько сгладить эф фект от грубости оценки максимального собственного числа.

Заметим, что при использовании (2.43) допускается неустойчи вость численной схемы на нескольких шагах. Однако это не приводит к катастрофическим последствиям, потому что при нарушении тре буемой точности за счет неустойчивости величина шага будет умень шена при проверке первого неравенства (2.42). Более того, интегриро вание с шагом, вызывающим неустойчивость численной схемы, может иметь и положительный эффект. Например, иногда возникает следую щая ситуация.

При интегрировании жестких задач явными методами на участке установления фактическая точность существенно превосходит зада ваемую точность. Как видно из формулы (2.35), это является следстви ем того, что на данном промежутке старые ошибки убывают за счет устойчивости численной схемы, а новые невелики в силу малости про изводных решения. В результате создается определенный запас точно сти, который можно использовать для повышения эффективности ал горитма интегрирования. Шаг можно увеличить до некоторого значения, не приводящего к нарушению (2.42), причем величина шага иногда в десятки раз превосходит максимально допустимую. После того как запас точности исчерпывается, шаг сбрасывается за счет кон троля точности до такого размера, что численная схема становится ус тойчивой.

При использовании (2.48) описанная выше ситуация достаточно часто реализуется в практических расчетах следующим образом. После прохождения переходного участка, соответствующего максимальному собственному числу матрицы Якоби, шаг устанавливается вблизи гра ницы области устойчивости hn D / | n max |, и начинает накапливать ся запас точности. Во многих случаях через некоторое время ошибка становится настолько мала, что по формуле (2.38) вместо значения n1 = n max вычисляется оценка n 2, где | n1 | | n 2 | | nN |, ni, 1 i N, – собственные числа матрицы Якоби системы (2.1).

В результате контроль устойчивости фактически прекращается. Так как переход от оценки n1 к n 2 осуществляется достаточно плавно – за несколько шагов, то и шаг интегрирования также увеличивается 2.4. Реализация явных методов плавно. Однако шаг может существенно вырасти, причем в силу фор мулы (2.43) он не будет уменьшен сразу после возникновения неус тойчивости. Уменьшение произойдет после того, как будет нарушено первое неравенство (2.42). В результате несколько шагов интегрирова ния будут выполнены с шагом, превышающим допустимый шаг. В та кой «взрывоопасной» ситуации имеет смысл контролировать оба нера венства (2.42), тем более что сама эта ситуация легко фиксируется.

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

Теперь сформулируем алгоритм интегрирования на основе явного метода с контролем точности и устойчивости численной схемы.

Шаг 1. Вычисляется оценка ошибки, p.

n Шаг 2. Вычисляется значение параметра q1 по формуле q1p ||, p || =.

n Шаг 3. Если q1 1, то величина шага hn полагается равной q1hn / 1,1 и происходит возврат на шаг 1.

Шаг 4. Вычисляется оценка ошибки, p.

n Шаг 5. Вычисляется значение параметра q2 по формуле p q2 ||, p || =.

n Шаг 6. Вычисляется оценка vn максимального собственного числа по формуле (2.38), причем если при всех j, 1 j N, выполняется j неравенство | n, p 1 |, где связана с разрядностью ЭВМ, то vn полагается равным 1.

Замечание. Случай vn = 1 означает, что оценка n max по форму ле (2.38) может вычисляться неправильно, поэтому требуется контро лировать оба неравенства (2.42).

Шаг 7. Если r 1 и q2 1, то величина шага hn полагается равной q2 hn / 1,1 и происходит возврат на шаг 1.

Г л а в а 2. КОНТРОЛЬ ТОЧНОСТИ И УСТОЙЧИВОСТИ ОДНОШАГОВЫХ МЕТОДОВ Шаг 8. Вычисляется приближение к решению по формуле (2.5) с шагом hn.

Шаг 9. Если q2 1, то величина прогнозируемого шага hn +1 пола гается равной hnew1, где hnew1 вычислен по формуле hnew1 = min(q1, q2 )hn. В противном случае новый шаг вычисляется по формуле (2.43).

Шаг 10. Выполняется следующий шаг интегрирования.

Описанный выше вычислительный процесс не следует рассматри вать как единственно возможный. В каждом конкретном случае реали зации могут различаться. Это может быть связано как со свойствами задач, на которые ориентируется алгоритм, так и с используемой чис ленной схемой. Поэтому ниже для некоторых численных формул алго ритм будет сформулирован отдельно.

Глава АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ С КОНТРОЛЕМ ТОЧНОСТИ НА ОСНОВЕ ЯВНЫХ МЕТОДОВ П остроены алгоритмы переменного шага на основе явных ме тодов с первого по четвертый порядок. Для контроля точно сти вычислений и при выборе величины шага интегрирования приме няется оценка ошибки (2.30).

3.1. ЯВНЫЙ МЕТОД ЭЙЛЕРА Во многих случаях расчеты требуется проводить с невысокой точ ностью – порядка 1% и ниже. Это связано с тем, что измерение кон стант, входящих в правую часть системы дифференциальных уравне ний, часто проводится достаточно грубо. Иногда такая точность расчетов удовлетворительна с точки зрения поставленной цели. Из вестно, что порядок аппроксимации численной схемы следует сочетать с требуемой точностью расчетов. Поэтому ниже будет рассмотрен ме тод первого порядка точности. В настоящее время многие пользовате ли применяют явный метод Эйлера с постоянным шагом для решения нежестких и умеренно жестких задач. Это обусловлено простотой реа лизации данной численной схемы, а также ее монотонностью. Отме тим, что для многих задач расчеты с постоянным шагом приводят к низкой эффективности алгоритма интегрирования. Ниже дан алгоритм переменного шага на основе явного метода Эйлера. Для решения (2.2) применим схему вида yn +1 = yn + hn f ( yn ), где hn – шаг интегрирова ния. Рассмотрение автономной задачи (2.2) связано с сокращением вы кладок, хотя построенный ниже алгоритм можно применять для неав тономных задач вида (2.1). В этом случае численная формула имеет вид yn+1 = yn + hn f (tn, yn ). Для построения алгоритма с контролем точности вычислений требуется оценить локальную ошибку и сформу Г л а в а 3. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ С КОНТРОЛЕМ ТОЧНОСТИ лировать алгоритм управления величиной шага. Далее предполагается, что алгоритм изменения величины шага построен на основе контроля точности численной схемы, а в неравенстве применяется оценка ло кальной ошибки n. Введем обозначения f n = f ( yn ), f n +1 = f ( yn +1 ) и построим неравенство для контроля точности.

Теорема 3.1. Пусть для численного решения задачи (2.2) использу ется явный метод Эйлера и пусть – требуемая точность расчетов.

Тогда для контроля точности вычислений имеет место неравенство 0,5h || f n +1 f n ||, где h = hn ;

f n и f n +1 – значения правой части системы (2.1) или (2.2) соответственно в точках tn и tn +1 ;

|| || – не которая норма в R N.

Используя разложения точного и приближенного решений в ряды Тейлора по степеням h до членов с h 2 включительно, получим, что для метода Эйлера выражение локальной погрешности n имеет вид n = 0,5h 2 f + O(h3 ), где элементарный дифференциал f f вычислен f на точном решении y (tn ). Учитывая разложение hf ( yn+1 ) в ряд Тейлора в окрестности приближенного решении yn до членов с h включительно hf ( yn +1 ) = hf n + h 2 f n f n + O(h3 ), получим h( f n +1 f n ) = = h 2 f n f n + O( h3 ). Сравнивая данное соотношение с локальной ошибкой, имеем неравенство для контроля точности вида 0,5h || f n +1 f n ||. Теорема доказана.

Отметим, что в случае успешного шага использование значения f n +1 в неравенстве для контроля точности к увеличению вычисли тельных затрат не приводит, потому что оно применяется на следую щем шаге для вычисления решения по схеме Эйлера. Теперь сформу лируем алгоритм интегрирования с контролем точности вычислений.

Пусть для численного решения задачи (2.1) используется явный метод Эйлера и пусть известно приближенное решение yn в точке tn с ша гом hn. Тогда для получения приближенного решения yn+1 в точке tn+1 справедлив следующий алгоритм.

Шаг 1. Вычислить приближение yn+1 по формуле yn+1 = = yn + hn f (tn, yn ).

3.2. Метод трапеций f n +1 = Шаг 2. Вычислить правую часть задачи (2.1), т. е.

= f (tn +1, yn +1 ).

|| n || = Шаг 3. Вычислить норму ошибки по формуле = 0,5h || f n +1 f n ||.

Шаг 4. Учитывая, что n = O(h 2 ), вычислить q по формуле q 2 || n || =.

Шаг 5. Если q 1, т. е. неравенство для контроля точности не вы полняется, перевычислить hn по формуле hn : = qhn / 1,1 и перейти на шаг 1.

Замечание. Постоянное число 1,1 введено в формулу для сокраще ния возвратов, т. е. повторных вычислений решения.

Шаг 6. Если q 1, т. е. вычисленное решение yn +1 в точке tn+ удовлетворяет требуемой точности, значение tn+1 принимается в каче стве следующего узла.

Шаг 7. Задать шаг hn +1 по формуле hn +1 = qhn / 1,1 для выполнения следующего шага интегрирования.

Шаг 8. Перейти на шаг 1.

Численные эксперименты на общепринятых тестовых задачах по казывают значительное повышение эффективности описанного алго ритма переменного шага по сравнению с расчетами с постоянной ве личиной шага.

3.2. МЕТОД ТРАПЕЦИЙ Для численного решения задачи Коши (2.1) рассмотрим метод вида zn +1 = yn + hn f (tn, zn ), (3.1) yn+1 = yn + 0,5hn [ f (tn, zn ) + f (tn +1, zn +1 ) ], z0 = y0.

Первая формула (3.1) – явный метод Эйлера, используемый для прогноза, вторая – формула трапеций, используемая для коррекции.

Заметим, что в формуле трапеций не используются итерации, а про гноз осуществляется с локальным порядком на единицу меньшим, чем коррекция. Так как имеется одновременно две формулы для вычисле ния решения, то контроль точности можно осуществлять без дополни тельных вычислений правой части системы (2.1). Будем интерпретиро Г л а в а 3. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ С КОНТРОЛЕМ ТОЧНОСТИ вать (3.1) как метод решения системы дифференциальных уравнений следующего вида:

z = f (t, z ), y = f (t, z ), t 0, z (0) = y (0) = y0. (3.2) Очевидно, что z (t ) y (t ).

Теорема 3.2. Пусть на решении z (t ) y (t ) задачи (3.2) выполнены условия | df / dt | A, | d 2 f / dt 2 | C, а в некоторой окрестно сти решения имеет место | f / z | B. Тогда zn и yn сходятся к решению y (t ) задачи (2.1) как O (ht2 ), где ht = max hk, 0 k n 1.

Для доказательства проинтегрируем уравнения (3.2) по переменной t на промежутке [tn, tn +1 ]. Вычитая из полученных выражений (3.1), имеем tn + z (tn +1 ) zn+1 = y (tn ) yn + [ f (t, z (t )) f (tn, zn )]dt, (3.3) tn y (tn +1 ) yn +1 = y (tn ) yn + tn + f (tn, zn ) + f (tn +1, zn +1 ) + f (t, z (t )) dt. (3.4) tn n = || z (tn ) zn || n = || y (tn ) yn || Введем обозначения и n hk, т. е.

|| z (tn ) yn ||. Здесь и далее будем предполагать, что tn = k = t0 = 0, tn+1 = tn + hn. Представим подынтегральные выражения в (3.3) и (3.4) в следующем виде:

f (t, z (t )) f (tn, zn ) = f ( t, z (t ) ) f (tn, z (tn )) + + f ( tn, z (tn ) ) f (tn, zn ), (3.5) f (tn, zn ) + f (tn +1, zn +1 ) f (t, z (t )) = f ( tn, z (tn ) ) + f ( tn +1, z (tn+1 ) ) = f (t, z (t )) f ( tn, z (tn ) ) f (tn, zn ) f ( tn+1, z (tn+1 ) ) f (tn +1, zn +1 ). (3.6) 2 3.2. Метод трапеций Используя в (3.5) и (3.6) формулу Тейлора с остаточным членом в форме Лагранжа первого и второго порядков соответственно, получим df f + ( z (tn ) zn ), f ( t, z (t ) ) f (tn, zn ) = (t tn ) dt z d2 f f (tn, zn ) + f (tn +1, zn +1 ) df f ( t, z (t ) ) = (t tn ) + (t tn )2 dt 2 dt 1 df 1 2 d 2 f f f + ( z (tn ) zn ) + ( z (tn +1 ) zn +1 ), hn hn z z 2 dt 4 dt где производные вычислены в соответствующих точках. Подставляя полученные выражения в (3.3) и (3.4) и переходя к норме, будем иметь n +1 n + Bhn n + Ahn / 2, (3.7) n +1 n + Bhn (n + n +1 ) / 2 + 5Chn / 12. (3.8) Введем обозначение rn = n + Bhn n. (3.9) Умножим (3.7) на выражение B (hn +1 + hn / 2) и сложим с (3.8). Усили вая полученное неравенство, имеем 3 rn +1 [1 + B (hn+1 + hn / 2)]rn + 5Chn / 12 + AB (hn +1 + hn / 2)hn / 2.

Учитывая, что при x 0 имеет место 1 1 + x e x, с использованием метода математической индукции легко показать следующее:

rn e B (tn +1 + 0,5tn ) ( D0tn + D1tn+1 ) ht2, (3.10) где D0 = 5C / 12 + AB / 4, D1 = AB / 2, ht = max hk.

0tn t, 0k n Из (3.9) и (3.10) следует, что n = || y (tn ) yn || = O(ht2 ). Покажем, что величина zn также сходится к y (tn ) с порядком O (ht2 ). Из (3.3) следует n +1 n + 2 Dhn, (3.11) Г л а в а 3. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ С КОНТРОЛЕМ ТОЧНОСТИ где D – максимум || f (t, z ) || в некоторой окрестности решения. Вычи тая первое уравнение (3.1) из второго, получим yn+1 zn+1 = hn [ f (tn+1, zn +1 ) f (tn, zn ) ] 2 = = hn f ( tn +1, z (tn +1 ) ) f ( tn, z (tn ) ) 2 + + hn f (tn +1, zn +1 ) f ( tn +1, z (tn+1 ) ) 2 hn f (tn, zn ) f ( tn, z (tn ) ) 2.

Используя в данном выражении формулу Тейлора с остаточным чле ном в форме Лагранжа первого порядка и переходя к норме, имеем || zn +1 yn+1 || Ahn / 2 + Bhn (n + n +1 ) / 2. (3.12) Подставим (3.11) в правую часть (3.12), т. е.

|| zn +1 yn +1 || Bhn ( n + n 1 ) / 2 + + Ahn / 2 + BDhn (hn + hn1 ). (3.13) Учитывая z (tn +1 ) zn+1 = z (tn +1 ) yn +1 + yn +1 zn +1, запишем n+ n+1 + | zn +1 yn+1 |. Так как n = O(ht2 ), то из этого неравенства и (3.13) следует n = O(ht2 ). Теорема доказана.

Теперь рассмотрим некоторые свойства численного решения, по лученного в результате применения схемы (3.1) к скалярному уравне нию y = y, Re( ) 0. В этом случае (3.1) приобретает следующий вид:

zn + zn =1 1.

y 1 + yn n +1 (1 + ) 2 Собственные числа матрицы перехода можно найти из уравнения 2 (3 / 2 + 1) + / 2 = 0, (3.14) а область устойчивости схемы (3.1) описывается кривой | | = 1 в ком плексной плоскости, Re( ) 0. В обозначениях = + i и 3.2. Метод трапеций = s + ig задача о нахождении области устойчивости сводится к на хождению множества точек, ограниченных кривой, 2 + 2 = 1. (3.15) Подставляя выражения для и в уравнение (3.14) и приравнивая действительную и мнимую части нулю, получим нелинейную систему алгебраических уравнений 2 2 (3s / 2 + 1) + 3 g / 2 + s / 2 = 0, 2 (3s / 2 + 1) 3 g / 2 + g / 2 = 0.

Учитывая соотношения (3.15), умножим первое уравнение на, а вто рое на и сложим. Затем умножим первое уравнение на, а второе на и сложим. Получим систему следующего вида:

(3s / 2 + 1) + ( s + gv) / 2 = 0, + 3 g / 2 + ( s g ) / 2 = 0.

Из второго уравнения выразим значение g = (2 s ) / (3 ), а затем, подставляя его в первое уравнение, получим s = s (), т. е.

s = 2( 1) 2 / (3 5), s (1) = 0, s (1) = 1. (3.16), т. е. имеем s Вычислим производную функции по s = ( 1)(6 14) / (3 5). Из полученной формулы видно, что на промежутке 1 1 функция s () монотонная, своего максимума достигает при = 1, т. е. 1 s 0. Выражая из (3.16), имеем = 1 + 3s / 4 ± 5,0625s 2 s. Так как s (1) = 1, то корень берем со знаком «минус». Подставляя в формулу для вычисления g при ус ловии (3.15), получим уравнение кривой ± 1 1 + 0,75s 0,5625s 2 s (2 s ) g=, { } 2 0,75s + 0,5625s s ограничивающей область устойчивости схемы (3.1) в плоскости.

Г л а в а 3. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ С КОНТРОЛЕМ ТОЧНОСТИ Для определения локальной ошибки выразим из уравнения (3.14) собственные числа матрицы перехода 13 1 + ± 0,5625 2 + + = 1,2 = 24 4 1 1 13 1 + ± + + 2 3 + O ( 4 ).

= 2 4 24 2 Разложим вектор ( zn, yn ) по собственным векторам данной матрицы:

1 / (1 ) zn = c1 + c2 = 2 yn 1 / (1 ) 1 c =. (3.17) 2 c Тогда вектор (c1, c2 ) можно представить в виде 2 zn c1 1 = 1.

2 1 c y 1 n Учитывая, что zn = yn + O (h 2 ), получим c1 = (1 )( 2 1) yn / ( 2 1 ) + O ( 2 ) = yn + O( ), (1 )[1 + 1 / (1 )] yn + O( 2 ) = c2 = 2 ( + 1 1 ) yn + O( 2 ) = O( 2 ).

= 2 Так как 2 = O ( ) и c2 = O ( 2 ), то из (3.17) следует, что вкладом в решение, даваемым вторым собственным вектором, можно пренеб речь. В результате относительную локальную ошибку схемы (3.1) в 3.3. Методы типа Рунге–Кутты точке tn +1 можно вычислить по следующей приближенной формуле:

n +1 = e 1 ( ) = 53 / 12 + O( 4 ). Тогда оценку глобальной ошибки схемы (3.1) можно вычислить по формуле (2.30), т. е. для контроля точности можно применять неравенство || yn +1 zn +1 ||.

3.3. МЕТОДЫ ТИПА РУНГЕ–КУТТЫ Ниже для численного решения задачи Коши (2.1) и (2.2) будем применять явные методы типа Рунге–Кутты следующего вида, соот ветственно i m yn+1 = yn + pmi ki, ki = hf tn + i h, yn + ij k j, (3.18) i =1 j = i m yn+1 = yn + pmi ki, ki = hf yn + ij k j, (3.19) i =1 j = где pmi, i, ij, 1 i m, 1 j i 1, – коэффициенты, определяю щие свойства точности и устойчивости (3.18) и (3.19), 1 = 0 ;

ki, 1 i m, – стадии метода. Там, где это не будет вызывать недоразуме ний, первый индекс при записи pmi, 1 i m, будем опускать. Для упрощения выкладок далее исследуем только методы (3.19). Однако все построенные ниже численные схемы можно использовать для ре шения неавтономных систем. Для этого достаточно в (3.18) подставить величины i, 1 i m, определенные следующим образом:

i i = j =1ij, 2 i m. (3.20) Запишем приближенное решение (3.19) в виде ряда Тейлора в ок рестности точки tn до членов с h4 включительно. Для этого введем в рассмотрение верхнюю треугольную матрицу Bm с элементами bij, которые имеют вид Г л а в а 3. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ С КОНТРОЛЕМ ТОЧНОСТИ b1i = 1, 1 i m, bki = 0, 2 k m, 1 i k 1;

(3.21) i bki = j = k 1ij bk 1, j, 2 k m, k i m, где ij есть коэффициенты схемы (3.19). В дальнейшем матрицы Bm будут использоваться для построения конкретных методов.

Теорема 3.3. Матрица Bm невырожденная тогда и только тогда, когда j, j 1 0, 2 j m, (3.22) где i, j 1 есть коэффициенты схемы (3.19).

Матрица Bm является верхней треугольной. С помощью (3.21) и с использованием метода математической индукции тривиально показы вается, что ее диагональные элементы можно представить в виде i bii = j =2 j, j 1, 2 i m. Доказательство теоремы следует из того, что условия bii 0, 2 i m, и (3.22) эквивалентны. Теорема до казана.

Соотношения (3.22) означают, что сразу после вычисления каждой новой стадии она вовлекается в вычислительный процесс, что является естественным условием при определении методов (3.18) или (3.19).

Поэтому всюду ниже будем считать условия (3.22) выполненными.

Теорема 3.4. Пусть метод (3.19) применяется для решения задачи = y, y (0) = y0, Re() 0. Тогда имеют место соотношения y m yn+1 = Qm ( z ) yn, z = h, Qm ( z ) = 1 + i =1 cmi z i, (3.23) m cmi = j =i bij pmj, 1 i m. (3.24) Сначала покажем справедливость равенства ( lj=1blj zl ) yn, 1 j m.

kj = (3.25) Доказательство проведем методом математической индукции. При j = 1 равенство (3.25) очевидно. Пусть (3.25) выполняется при j = i.

Тогда для j = i + 1 имеем 3.3. Методы типа Рунге–Кутты )( ) ( i i j ki +1 = z yn + j =1i +1, j k j = z 1 + j =1 i +1, j l =1blj z l yn.

Учитывая b1, i +1 = 1 и уравнения (3.21), запишем выражение в скоб ках по возрастающим степеням z, т. е. ki +1 = b1,i +1 z + l =1 z l + i ) y. Отсюда, с учетом последнего равенства (3.21), ( i b j =l i +1, j lj n ( ) ki +1 = b1,i +1 z + l =1bl +1,i +1 z l +1 yn. Теперь, заменяя индекс l + 1 на l и i внося первое слагаемое под знак суммы, получим (3.25). Для доказа тельства соотношения (3.23) подставим (3.25) в первую формулу ( ) (3.19). В результате имеем yn+1 = 1 + i =1 pmi l =1bli z l yn. Пере m i писывая выражение в скобках по степеням значения z и учитывая ра венства (3.21), получим (3.23), (3.24). Теорема доказана.


Функция устойчивости Qm ( z ) явного m -стадийного метода типа Рунге–Кутты представляет собой многочлен степени m относительно z. Данная теорема позволяет вычислить коэффициенты этого полино ма через коэффициенты (3.19). Поэтому перепишем их еще в несколь ко ином виде. Введем обозначения Cm = (cm1,…, cmm )T и Pm = ( pm1,…, pmm )T. Тогда (3.24) можно представить в виде Bm Pm = Cm, (3.26) где элементы матрицы Bm определены формулами (3.21).

В случае задачи Коши для линейной системы обыкновенных диф ференциальных уравнений с постоянными коэффициентами y = Ay + b, y (t0 ) = y0, t0 t tk, (3.27) соотношения (3.26) можно применять для построения методов задан ного порядка точности. Используя обозначение f ( y ) = Ay + b, нетруд но убедиться, что для точного решения y (tn +1 ) задачи (3.27) и при ближенного решения yn +1, полученного по схеме (3.19), справедливы представления y (tn+1 ) = y (tn ) + i =1 hi f (i 1) f / i!, (3.28) Г л а в а 3. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ С КОНТРОЛЕМ ТОЧНОСТИ yn+1 = yn + i =1 cmi hi f n(i 1) f n, m (3.29) где элементарные дифференциалы вычислены соответственно на точ ном y (tn ) и приближенном yn решениях, а коэффициенты cmi, 1 i m, и коэффициенты схемы (3.19) связаны равенствами (3.26).

Из сравнения (3.28) и (3.29) при условии yn = y (tn ) следует, что p -й порядок точности, если схема (3.19) будет иметь cmi = 1 / i!, 1 i p. Теперь, задаваясь, например, значениями парамет ров ij, 2 i m, 1 j i 1, остальные коэффициенты схемы (3.19) p -го порядка точности получим из линейной системы (3.26). Таким образом, если метод (3.19) применяется для решения задачи (3.27), то справедливы следующие следствия из теоремы 3.4.

Следствие 1. Для любого целого положительного m можно по строить m -стадийный метод m -го порядка точности.

Следствие 2. Максимальный порядок точности m -стадийной схемы не может быть больше m.

Запишем представление точного решения y (tn +1 ) в окрестности точ ки tn в виде ряда Тейлора по степеням h до членов с h4 включительно.

Введем следующие обозначения элементарных дифференциалов:

N N f N fi f f = f 2 f = fi, f j, f i =1 yi i =1 yi j =1 y j f N fi N f j 2 f NN N f 2 = f 3 f = fi f j, fk, f i =1 j =1 yi y j i =1 yi j =1 y j k =1 yk f N N 2 fi N f f 2 = f f j fk, (3.30) i =1 yi j =1 k =1 yi yk 2 f N f NN f f 2 = i fk f j, f i =1 j =1 yi y j k =1 yk 3 f NN N f 3 = y y fi f j f k.

f i j yk i =1 j =1 k = 3.3. Методы типа Рунге–Кутты Тогда точное решение y (tn +1 ) задачи (2.2) в окрестности точки tn можно записать следующим образом:

y (tn +1 ) = y (tn ) + hf + h 2 f / 2 + h3 ( f 2 f + f 2 ) / 6 + f f (3.31) + h 4 ( f 3 f + f f 2 + 3 f 2 + f 3 ) / 24 + O(h5 ), f ff f где элементарные дифференциалы вычислены на точном реше нии y (tn ).

Приближенное решение yn+1, полученное по (3.19), представимо в виде m m m yn+1 = yn + bij pmj hi f n(i 1) f n + b2i pmi h3 f n f n2 2 + i =2 i =1 j =i m 3 m + h 4 b2i pmi f n f n + 6 b2i b3i pmi f n f n f n2 + 3 i = 2 i =3 m i 1 2 + 6 ij b2 j pmi f n f n f n2 6 + O( h5 ), (3.32) i =3 j = 2 где индекс n при записи элементарных дифференциалов (3.30) означа ет, что они вычислены в точке yn, а коэффициенты bij, 1 i, j m, определены формулами (3.21). Справедливость (3.32) вытекает из сле дующей теоремы.

Теорема 3.5. Стадии метода (3.19) представимы в виде i b ji h j fn( j 1) fn + b2i h3 fn fn2 / 2 + h4 b2i f nf n3 + 2 ki = j = i 1 2 + 6b2i b3i f n f n f n2 + 3 ij b2 j f n f n f n2 6 + O (h5 ).

(3.33) j =2 Правильность представления (3.33) тривиально доказывается мето дом математической индукции. Из-за громоздкости выкладок доказа тельство не приводится. Подставляя (3.33) в первую формулу (3.19) и учитывая верхний треугольный вид матрицы Bm, получим (3.32).

Г л а в а 3. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ С КОНТРОЛЕМ ТОЧНОСТИ 3.4. МЕТОДЫ ТИПА РУНГЕ–КУТТЫ ВТОРОГО ПОРЯДКА ТОЧНОСТИ Для численного решения задачи Коши (2.2) рассмотрим двухста дийную формулу типа Рунге-Кутты yn+1 = yn + p1k1 + p2 k2, k1 = hf ( yn ), k2 = hf ( yn + 21k1 ). (3.34) Получим условия второго порядка точности схемы (3.34). Для этого запишем приближенное решение в виде ряда Тейлора. С использова нием (3.32) имеем yn+1 = yn + ( p1 + p2 )hf n + 2 p2 h 2 f n f n + 2 p2 h3 f n f n2 / 2 + O(h 4 ), где параметр 2 задается формулой (3.20). Пусть yn = y (tn ). Сравни вая полученное представление приближенного решения yn +1 и (3.31) до членов с h 2 включительно, получим условия второго порядка, т. е.

p1 + p2 = 1, 2 p2 = 1 / 2, (3.35) при этом локальная ошибка n,2 схемы (3.34) имеет вид ( ) n,2 = h3 f 2 f + 1 3 2 p2 f 2 6 + O (h 4 ).

f Рассмотрим одновременно метод первого порядка zn +1 = yn + k1, (3.36) локальная ошибка которого имеет вид n,1 = h 2 f / 2 + O (h3 ). Условие f согласования (2.22) будет выполнено, если 2 p2 = 1 / 3. Теперь коэф фициенты схемы (3.34) однозначно определяются из (3.35) и имеют вид 2 = 21 = 2 / 3, p1 = 1 / 4, p2 = 3 / 4. (3.37) Введем обозначение n,2 = ( yn zn ) / 6 (k2 k1 ) / 4 = h 2 f n f n / 6 + +O (h3 ). Тогда для контроля точности вычислений можно использо вать неравенство типа (2.31) || k2 k1 || 4. (3.38) 3.4. Методы типа Рунге–Кутты второго порядка точности Отметим, что приращение k1 зависит от размера шага линейно, поэто му повторное вычисление решения в случае нарушения требуемой точности расчетов – при возврате в результате невыполнения неравен ства (3.38) – будет приводить только к одному дополнительному вы числению правой части исходной задачи (2.2). Учитывая, что hf ( yn +1 ) k1 = h 2 f n f n + O (h3 ), при выборе шага будем контролировать оценку ошибки вида n,2 = [hf ( yn+1 ) k1 ] / 6, т. е. при выборе шага будем дополнительно проверять неравенство || hf ( yn +1 k1 || 6. (3.39) На рис. 3.1–3.5 слева приведены линии уровня | Q3,2 ( z ) |= s при значениях s, равных числам 1, 0,7;

0,3;

0,1;

справа – объемные изо бражения областей устойчивости. При сравнении свойств устойчиво сти численных схем удобно пользоваться отношением длины интерва ла устойчивости к количеству вычислений функции f на шаге интегрирования. В явном методе Эйлера, например, на одно вычисле ние функции f приходится две единицы длины интервала устойчиво сти. Исследуем схему (3.34) на линейной задаче y = y. Применяя (3.34), (3.37) для решения данной задачи, получим yn +1 = Q2,2 ( z ) yn, Q2,2 ( z ) = 1 + z + z 2 / 2. Схема (3.34) будет устойчива, если z = h, | Q2,2 ( z ) | 1. Если переменная z вещественная, то это неравенство равносильно 1 + 0,5 z 0, z 2 + 2 z + 4 0, т. е. 2 z 0. Отсюда видно, что в алгоритме RK21 на одно вычисление функции f приходится единица длины интервала устойчивости. Если шаг ограничен только в силу устойчивости, то явный метод Эйлера может оказаться эффек тивнее. Область устойчивости формулы (3.34) приведена на рис. 3.1.

Построим еще одну численную схему второго порядка с контролем точности вычислений, в которой на одно вычисление правой части ис ходной задачи (2.2) приходится не менее двух единиц длины интервала устойчивости. Для этого рассмотрим формулу типа Рунге–Кутты с тремя вычислениями функции f на шаге интегрирования yn +1 = yn + p1k1 + p2 k2 + p3k3, k1 = hf ( yn ), (3.40) k2 = hf ( yn + 21k1 ), k3 = hf ( yn + 31k1 + 32 k2 ).

Г л а в а 3. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ С КОНТРОЛЕМ ТОЧНОСТИ Рис. 3.1. Область устойчивости метода RK Используя (3.32), получим следующее представление приближен ного решения yn+1 = yn + ( p1 + p2 + p3 )hf n + ( 2 p2 + 3 p3 )h 2 f n f n + ( ) + h3 gf n2 f n + 2 p2 + 3 p3 f n f n2 / 2 6 + O( h3 ), 2 (3.41) где g = 232 p3, а i, 2 i 3, задаются по формуле (3.20).

Пусть yn = y (tn ). Сравнивая (3.31) и (3.41) до членов с h2 включи тельно, получим условия второго порядка точности схемы (3.40), т. е.

p1 + p2 + p3 = 1, 2 p2 + 3 p3 = 1 / 2. (3.42) Используя для контроля точности формулу (3.36), видим, что для вы полнения условий согласования (2.22) достаточно потребовать выпол нения соотношения 2 p2 + 3 p3 = 1 / 3, (3.43) при этом локальная ошибка n,2 формулы (3.40) имеет вид n,2 = (1 6 g )h3 f 2 f / 6 + O(h 4 ). (3.44) Из (3.44) следует, что при выполнении (3.42), (3.43) и равенства g = 1 / 6 схема (3.40) имеет третий порядок. Выбирая g близким к 1 / 6, можно коэффициент в главном члене локальной ошибки сделать каким угодно малым. Однако это может приводить к доминированию слагаемых порядка O (h 4 ), и характер поведения ошибки уже не будет определяться главным членом.

3.4. Методы типа Рунге–Кутты второго порядка точности Исследуем совместность (3.42), (3.43) при условии g = 232 p3.

Сначала рассмотрим следующую систему:

2 2 p2 + 3 p3 = 1 / 2, 2 p2 + 3 p3 = 1 / 3, 232 p3 = g, линейную относительно p2 и p3. Для ее совместности необходимо обращение в нуль определителя расширенной матрицы. Из (3.20) и (3.22) следует, что 2 0. Раскрывая определитель и приводя подоб ные члены, получим 6 g 3 (3 2 ) = 32 2 (2 3 2 ). (3.45) Выберем g, 2, 3 и 32 так, чтобы они удовлетворяли (3.45), причем положим 2 = 21, 3 = 31 + 32. Тогда коэффициенты pi, 1 i 3, можно найти последовательным исключением из линейной системы:

232 p3 = g, 2 p2 + 3 p3 = 1 / 2, p1 + p2 + p3 = 1. (3.46) Учитывая (3.36), для контроля точности вычислений (3.40) можно ис пользовать оценку ошибки n,2 вида n,2 = (1 6 g )[( p1 1)k1 + p2 k2 + p3k3 ] / 3 = = (1 6 g )h 2 f n f n / 6 + O( h3 ), (3.47) т. е. на каждом шаге нужно контролировать неравенство || ( p1 1)k1 + p2 k2 + p3k3 || 3 / | 1 6 g |. (3.48) При применении (3.48) в случае невыполнения точности расче тов повторное вычисление решения будет сопровождаться двумя дополнительными вычислениями правой части исходной задачи (2.1). Поэтому построим более экономичное неравенство для кон троля точности вычислений. Сравнивая (3.47) с соотношением k2 k1 = 2 h 2 f n f n + O( h3 ), видим, что для контроля точности (3.40) с коэффициентами (3.45), (3.46) можно использовать неравенство || k2 k1 || 6 | 2 | / | 1 6 g |. (3.49) Так как k1 зависит от размера шага линейно, то в этом случае возврат при нарушении требуемой точности расчетов приводит всего лишь к одному дополнительному вычислению функции f.


Г л а в а 3. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ С КОНТРОЛЕМ ТОЧНОСТИ Дополнительное неравенство для выбора величины шага интегри рования построим так же, как и для схемы (3.34). Учитывая вид разно сти hf ( yn +1 ) k1, оценку ошибки 2 вычислим по формуле n,,2 = | 1 6 g | [hf ( yn +1 ) k1 ] / 6, т. е. при выборе шага, наряду с (3.49), n будем контролировать неравенство h || f ( yn+1 ) f ( yn ) || 6 / | 1 6 g |. (3.50) Приведем некоторые соображения по выбору конкретных коэффи циентов схемы (3.40). Из (3.45) и (3.46) видно, что при их выборе име ется некоторый произвол. Свободные параметры g, 2, 3 и можно использовать, например, для улучшения свойств устойчивости схемы (3.40). Для определения влияния величины g на размер интер вала устойчивости применим (3.40) при решении тестового уравнения y = y. Учитывая (3.45) и (3.46) с применением (3.23) и (3.24), полу чим функцию устойчивости Q3,2 ( z ) схемы (3.40), т. е.

Q3,2 ( z ) = 1 + z + z 2 / 2 + gz 3, z = h.

Вопрос о том, как связаны коэффициенты многочлена устойчиво сти (3.23) с размером и формой области устойчивости, ниже будет рас сматриваться подробно. Здесь ограничимся замечанием, что макси мальная длина 6, 26 интервала устойчивости достигается при g = 1 / 16. Однако при данном значении g имеет место Q3,2 (4) = 1, и поэтому небольшие возмущения, например, за счет ошибок округле ний, могут приводить к уменьшению области устойчивости. Если ис ходить из требования, чтобы длина интервала устойчивости была дос таточно велика, то значение g можно выбрать равным 1 / 15. При этом на одно вычисление правой части задачи (2.1) приходится примерно две единицы длины интервала устойчивости. С точки зрения устойчи вости данный метод и метод Эйлера различаются незначительно. В то же время формула (3.40) обладает вторым порядком точности и снаб жается алгоритмом контроля точности вычислений.

Если Q3,2 ( z ) обладает монотонностью на интервале устойчивости, а таким свойством обладает известный метод Мерсона, то можно ожи дать, что небольшие возмущения не приведут к сокращению области устойчивости. Требование монотонности эквивалентно тому, что 3.4. Методы типа Рунге–Кутты второго порядка точности производная Q3,2( z ) = 3gz 2 + z + 1 не меняет знак на интервале устой чивости. Учитывая, что корни уравнения Q3,2 ( z ) = 0 имеют вид z1,2 = g 1 6 ± g 1 (1 12 g ) 6, убеждаемся, что Q3,2 ( z ) будет моно тонной, если g 1 / 12.

Далее в процессе вычислений по формуле (3.40) приходится опре делять приближенное решение в двух промежуточных точках yn+11 = yn + 21k1 и yn+1,2 = yn + 31k1 + 32 k2. Если области устойчи, вости промежуточных численных схем же области устойчивости ос новной схемы, то приближение к решению в промежуточных точках может искажаться за счет усиления ошибок округлений, что приводит к дополнительной ошибке в решении. Потребуем, чтобы интервалы устойчивости этих схем были не меньше, чем у основной схемы.

Нетрудно убедиться, что при 21 = 1 / 3 и 3 = 232 длины интерва лов устойчивости промежуточных численных схем не меньше шести.

Рис. 3.2. Область устойчивости метода (3.40), (3.51)      Рис. 3.3. Область устойчивости метода (3.40), (3.52) Г л а в а 3. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ С КОНТРОЛЕМ ТОЧНОСТИ Теперь, полагая значение параметра g равным 1 / 12, 1 / 15 и 1 / 16, из (3.45) и (3.46) получим следующие наборы коэффициентов:

1 g=, 2 = 21 = 31 = 32 =, 12 (3.51) 2 1 3 =, p1 =, p2 = 0, p3 = ;

3 4 1 1 3 g =, 2 = 21 =, 3 =, 31 = 32 =, 15 3 4 (3.52) 1 3 p1 =, p2 =, p3 = ;

6 10 1 g =, 2 = 21 =, 31 = 32 =, 16 3 7 1 3 3 =, p1 =, p2 =, p3 =.

9 7 8 Области устойчивости методов (3.40) с коэффициентами (3.51) и (3.52) приведены на рис. 3.2 и 3.3.

3.5. МЕТОДЫ ТИПА РУНГЕ–КУТТЫ ТРЕТЬЕГО ПОРЯДКА ТОЧНОСТИ Легко показать, что на основе трехстадийной схемы (3.40) нельзя построить метод третьего порядка с контролем ошибки (2.30). Поэтому рассмотрим четырехстадийную схему yn+1 = yn + pi ki, k1 = hf ( yn ), k2 = hf ( yn + 21k1 ), i = k3 = hf ( yn + 31k1 + 32 k2 ), (3.53) k4 = hf ( yn + 41k1 + 42 k2 + 43k3 ).

С использованием (3.32) имеем yn+1 = yn + pi hf n + ( 2 p2 + 3 p3 + 4 p4 )h 2 f n f n + i = h3 ( 232 p3 + 242 p4 + 343 p4 ) f n 2 f n + 3.5. Методы типа Рунге–Кутты третьего порядка точности ( ) + 2 p2 + 3 p3 + 4 p4 f n f n2 / 2 + 2 2 ( ) + h 4 23243 p4 f n3 f n + 232 p3 + 242 p4 + 343 p4 f n f n f n2 / 2 + 2 2 ( ) + 3 p2 + 3 p3 + 3 p4 f n f n / 6 + O( h5 ), 3 2 где элементарные дифференциалы вычислены на приближенном ре шении yn, а параметры i, 2 i 4, определены формулами (3.20).

Для контроля точности (3.53) будем использовать формулу второго порядка точности вида zn +1 = yn + b1k1 + b2 k2, (3.54) где стадии k1 и k2 заданы в (3.53).

Проводя рассуждения, аналогичные тем, которые применялись при исследовании (3.40), получим условия третьего порядка точности схе мы (3.53) и условия согласования ошибок (2.22), т. е.

1) p1 + p2 + p3 + p4 = 1;

2) 2 p2 + 3 p3 + 4 p4 = 1 / 2;

3) 2 p2 + 3 p3 + 2 p4 = 1 / 3;

4) 3 p2 + 3 p3 + 3 p4 = 1 / 4;

2 2 4 2 5) 232 p3 + 242 p4 + 343 p4 = 1 / 6;

(3.55) 6) 2332 p3 + 2 442 p4 + 3 443 p4 = 1 / 8;

7) 232 p3 + 242 p4 + 3 43 p4 = {(24 g 1)(2 3 2 ) + 2} / 24;

2 8) g = 23243 p4, где g – некоторое число, g 0. При записи седьмого соотношения (3.55) воспользовались равенством b2 2 = 1 / 2, которое следует из ус ловия второго порядка точности схемы (3.54).

Исследуем совместность (3.55). Если выбрать параметры g, i, ij, 2 i 4, 1 j i 1, то значения коэффициентов рi, 1 i 4, можно вычислить последовательным исключением из системы 23243 p4 = g, 232 p3 + 242 p4 + 343 p4 = 1 / 6, (3.56) 2 p2 + 3 p3 + 4 p4 = 1 / 2, p1 + p2 + p3 + p4 = 1.

Г л а в а 3. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ С КОНТРОЛЕМ ТОЧНОСТИ Возьмем пятое, седьмое и восьмое уравнения (3.55). Если их рассмат ривать как систему линейных уравнений относительно p3 и p4, то для ее совместности требуется обращение в нуль определителя расши ренной матрицы. После приведения подобных членов это условие име ет вид 24 g 3 (3 2 ) + 232 [4 2 (24 g 1)(2 3 2 ) 2] = 0. (3.57) Теперь умножим второе, третье и пятое уравнения (3.55) на 3 и вы чтем из полученных соотношений третье, четвертое и шестое равенст ва соответственно. Будем иметь 1) 4 / 2 1 / 3 = 2 ( 4 2 ) p2 + 3 ( 4 3 ) p3 ;

2) 4 / 4 1 / 4 = 2 ( 4 2 ) p2 + 3 ( 4 3 ) p3 ;

(3.58) 3) 4 / 6 1 / 8 = 232 ( 4 3 ) p3.

Умножим первое уравнение (3.58) на 2 и вычтем второе. В результа те имеем ( 4 / 3 1 / 4) 2 ( 4 / 2 1 / 3) = 3 (3 2 )( 4 3 ) p3.

Подставим сюда вместо ( 4 3 ) p3 его выражение из последнего со отношения (3.58). Окончательно получим 3 (3 2 )(4 4 3) = 232 (8 2 + 8 4 12 2 4 6). (3.59) Далее, из последнего соотношения (3.55) следует, что 2 и 32 от личны от нуля. Тогда 4 3 / 4, так как в противном случае из (3.59) получим, что либо 2 = 0, либо 32 = 0. Умножим (3.57) на соотно шение (4 4 3), а (3.59) – на 24g. Складывая полученные выраже ния, имеем 4 = 6 g + 3 / 4.

Потребуем совместности второго, третьего, четвертого и восьмого уравнений (3.55). Это будет выполнено, если равен нулю определитель соответствующей расширенной матрицы. Раскрывая определитель, получим 23243 (12 2 3 83 8 2 + 6) = 24 g 4 ( 4 3 )( 4 2 ).

3.5. Методы типа Рунге–Кутты третьего порядка точности Наконец потребуем совместности пятого, шестого и восьмого уравне ний (3.55). Это требование будет выполнено, если обращается в нуль определитель расширенной матрицы, что равносильно выполнению равенства 24 g 242 ( 4 3 ) = 23243 (3 43 ) + 24 g 343 (3 4 ).

Если выбрать параметры g, ij, 2 i 4, 1 j i 1, из системы уравнений 2 = 21;

3 = 31 + 32 ;

4 = 41 + 42 + 43 ;

4 = g + 1 / 8;

232 (8 2 + 8 4 12 2 4 6) = 24 g 3 (3 2 );

(3.60) 23243 (12 2 3 83 8 2 + 6) = 24 g 4 ( 4 3 )( 4 2 );

g 242 ( 4 2 ) = 43 [ 232 (3 43 ) + 24 g 3 (3 4 ) ] 24, то затем pi, 1 i 4, находятся последовательным исключением из (3.56). В результате получим метод третьего порядка с локальной ошибкой n,3 = (1 24 g )h 4 f f 2 f + (2 3 2 ) f 2 / 2 24 + O (h5 ). (3.61) f Учитывая, что локальная ошибка n,2 формулы (3.54) имеет вид n,2 = h3 f 2 f + (2 3 2 ) f 2 / 2 6 + O( h 4 ), f для контроля точности можно применять неравенство вида (2.31), т. е.

|| yn zn || 4 / | 1 24 g |, (3.62) где yn и zn определены соответственно формулами (3.53) и (3.54). Так как yn zn = ( p1 b1 ) k1 + ( p2 b2 ) k2 + p3k3 + p4 k4, то неравенство (3.62) можно переписать в виде || ( p1 b1 )k1 + ( p2 b2 )k2 + p3k3 + p4 k4 || 4 / | 1 24 g |, где pi, 1 i 4, определены в (3.56), а b1 и b2 вычисляются из требо вания второго порядка точности схемы (3.54), т. е. b2 = 0,5 / 2 и b1 = 1 b2.

Г л а в а 3. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ С КОНТРОЛЕМ ТОЧНОСТИ Теперь приведем некоторые соображения по выбору свободных коэффициентов в системе (3.60). Сначала потребуем, чтобы в главный член локальной ошибки входило только слагаемое вида h 4 f 3 f. Оче видно, что это условие эквивалентно равенству 2 = 2 / 3. Далее по ложим 3 = 1. В этом случае стадия k3 вычисляется в точке tn+1, и отпадает необходимость в построении дополнительного неравенства для выбора величины шага интегрирования. И наконец последний сво бодный параметр g используем для увеличения области устойчивости (3.53). Применяя (3.53) для численного решения тестовой задачи y = y, Re( ) 0, получим (3.23), где функция устойчивости Q4,3 ( z ) имеет вид Q4,3 ( z ) = 1 + z + z 2 / 2 + z 3 / 6 + gz 4, z = h. Ясно, что длина интервала устойчивости схемы (3.53) зависит от выбора параметра g.

С помощью ЭВМ установлено, что максимальная длина интервала устойчивости достигается при значении g 1 / 54 и равна примерно шести. Однако так как при g = 1 / 54 имеет место Q4,3 (4,5) 1, Рис. 3.4. Область устойчивости метода (3.53), (3.63) (показаны линии уровня | Q3,2 ( z ) | = s при s, равном 1;

0,7;

0,3;

0,1) Рис. 3.5. Область устойчивости метода (3.53), (3.64) 3.6. Метод Рунге–Кутты–Мерсона то небольшие возмущения могут приводить к уменьшению области устойчивости. Поэтому для повышения надежности расчетов имеет смысл выбрать g = 1 / 53.

Далее, при g = 1 / 48 коэффициент в локальной ошибке (3.61) меньше, чем при g = 1 / 53, а в некоторой окрестности точки z = 4, выполняется неравенство | Q4,3 ( z ) | 0,5. Такой запас может быть по лезен при небольших возмущениях, например, за счет ошибок округ лений.

Итак, пусть 2 = 2 / 3 и 3 = 1, а параметр g принимает значения 1 / 48 и 1 / 53. Тогда из (3.53) и (3.56) последовательным исключением найдем следующие два набора коэффициентов:

1 2 11 g = ;

2 = 21 = ;

3 = 1;

31 = ;

32 = ;

48 3 8 7 1351 525 4 = ;

41 = ;

42 = ;

43 = ;

(3.63) 8 1024 1024 17 27 2 p1 = ;

p2 = ;

p3 = ;

p4 = ;

84 20 3 1 2 71 18 g= ;

2 = 21 = ;

3 = 1;

31 = ;

32 = ;

4 = ;

53 3 53 53 24387129 9264375 41 = ;

42 = ;

43 = ;

(3.64) 19056256 19056256 443 693 53 p1 = ;

p2 = ;

p3 = ;

p4 =.

2196 500 87 Области устойчивости методов (3.53) с коэффициентами (3.63) и (3.64) приведены на рис. 3.4 и 3.5.

3.6. МЕТОД РУНГЕ–КУТТЫ–МЕРСОНА Одним из самых эффективных среди явных методов типа Рунге– Кутты четвертого порядка точности является метод Мерсона, который имеет вид 1 2 yn+1 = yn + k1 + k4 + k5 ;

k1 = hf ( yn );

6 3 Г л а в а 3. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ С КОНТРОЛЕМ ТОЧНОСТИ 1 k2 = hf yn + k1 ;

k3 = hf yn + k1 + k2 ;

(3.65) 3 3 1 1 k4 = hf yn + k1 + k3 ;

k5 = hf yn + k1 k3 + 2k4.

8 2 Пятое вычисление правой части не дает увеличения порядка точности до пятого, однако позволяет расширить интервал устойчивости при мерно до величины 3,5 по вещественной оси. Кроме того, что более существенно, пятое вычисление правой части задачи (2.1) или (2.2) позволяет оценить величину локальной ошибки n,4 с помощью ста дий ki, 1 i 5, т. е.

n,4 = (2k1 9k3 + 8k4 k5 ) / 30. (3.66) Во многих библиотеках программ имеется программная реализация метода Мерсона со следующим алгоритмом выбора величины шага интегрирования.

Шаг 1. Вычисляются приращения ki, 1 i 5.

Шаг 2. Вычисляется оценка локальной ошибки по формуле (3.66).

Шаг 3. Если для некоторого значения j, 1 j N, выполняется неравенство j j | n,4 | | yn |, (3.67) где j есть номер компоненты, то шаг уменьшается в 2 раза и снова вычисляются ki, 1 i 5, т. е. управление передается на шаг 1.

Замечание. Если для некоторого j, 1 j N, выполняется нера j венство | yn | Rzero, то по этой компоненте условие (3.67) не прове ряется. Величина Rzero определяется длиной слова на компьютере.

Шаг 4. Вычисляется решение по формуле (3.65).

Шаг 5. Если для всех j, 1 j N, выполняется условие j j | n,4 | | yn | /32, (3.68) то шаг удваивается. В противном случае остается без изменения.

Шаг 6. Выполняется следующий шаг интегрирования.

В формуле (3.68) используется константа 1 / 32, которая вводится следующим образом. Пусть после успешного шага попытаемся увели 3.6. Метод Рунге–Кутты–Мерсона чить шаг в 2 раза, т. е. hn +1 = 2hn. Учитывая, что n,4 = O( h5 ), запи шем n,4 ( hn +1 ) = n,4 (2hn ) = 25 n,4 (hn ), откуда и следует (3.68). Изме нение шага интегрирования в 2 раза обладает некоторыми недостатка ми. Если при вычислении решения шаг интегрирования часто и резко меняется, например, при большом количестве переходных процессов, такой алгоритм может приводить к неоправданно большим затратам времени. Кроме того, при неудачном выборе начального шага может понадобиться достаточно много шагов для достижения максимального шага.

Обоснование оценки ошибки для контроля точности вычислений проведем на линейном уравнении y = y, y (0) = y0, Re() 0. При меняя (3.65) для решения этой задачи, получим (3.23), где функция роста Q5,4 ( z ) имеет вид Q5,4 ( z ) = 1 + z + z 2 / 2 + z 3 / 6 + z 4 / 24 + z 5 / 144, z = h.

Учитывая, что exp( z ) = z i / i, получим, что относительная локальная i = n, ошибка численной формулы (3.65) имеет вид 5 n,4 = z / 720 + O( z ). Используя формулу (2.30), получим оценку ошибки n,4, т. е. n,4 z 4 / 720. Пусть заданная относительная точ ность вычислений равна. Тогда | z | 4 720. Теперь можно запи сать неравенство для контроля точности через оценку локальной ошибки, т. е. | n,4 | 55/ 4, при получении которого применялось со отношение 55/ 4 7201/ 4 5/ 4. Несмотря на то что обоснование данно го неравенства проведено на линейном скалярном уравнении, оно с достаточно высокой надежностью использовалось для решения нели нейных задач. Кроме того, неравенство такого типа успешно применя лось в алгоритмах интегрирования на основе других численных схем.

Заметим, что оценка (3.66) является величиной О(h5 ) только в случае линейной задачи, т. е. когда имеет место соотношение f (t, y ) = a + bt + cy. В других случаях n,4 содержит некоторые члены Г л а в а 3. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ С КОНТРОЛЕМ ТОЧНОСТИ порядка O(h 4 ). Разложением (3.65) в ряд Тейлора до членов с h включительно легко убедиться, что для скалярного уравнения (2.2) выполняется соотношение n,4 = h 4 [2 f 3 + 6 f f 2 + 3 f 2 ] / 45 + O(h5 ).

f f ff При построении алгоритма интегрирования этот факт можно учи тывать следующим образом. Предположим, что n,4 ведет себя при увеличении шага как O(h5 ), а при уменьшении как O(h 4 ). Такой небольшой запас в точности позволяет частично избежать возможных повторных вычислений решения. Так как k5 вычисляется в точке tn +1, то отпадает необходимость в построении дополнительного неравенст ва для выбора величины шага интегрирования.

Глава АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ С КОНТРОЛЕМ УСТОЙЧИВОСТИ ЧИСЛЕННОЙ СХЕМЫ В данной главе построены неравенства для контроля устойчи вости явных методов типа Рунге–Кутты со второго по пятый порядок точности. Ниже для контроля устойчивости явных методов типа Рунге–Кутты (3.19) будем использовать неравенство вида (2.41).

Для получения оценки максимального собственного числа требуется выполнение требования типа (2.37). Напомним, что параметр p в ус ловиях согласования локальных ошибок (2.37) определяет число ите раций в степенном методе. Поэтому требования типа (2.37) можно рас сматривать при некотором значении l, l p. За счет этого можно избежать дополнительных соотношений на параметры численной схе мы и тем самым избавиться от дополнительных вычислений правой части дифференциальной задачи. Поэтому для методов (3.19) посту пим следующим образом.

Выберем постоянные коэффициенты bi и bi, 1 i m + 1, из усло вия выполнения соотношений m bi ki + bm+1hf ( yn+1 ) = hl +1 f nl f n + O(hl +3 ), i = (4.1) m biki + bm+1hf ( yn+1 ) = hl + 2 f n(l +1) f n + O(hl +3 ), i = где ki, 1 i m, определены в (3.19), а элементарные дифференциалы вычислены в точке yn. Для получения (4.1) следует разложить f ( yn+1 ) в ряд Тейлора и подставить это представление и (3.33) в ле вую часть (4.1). Затем нужно приравнять нулю соотношения на пара метры bi и bi, 1 i m + 1, перед соответствующими слагаемыми.

Г л а в а 4. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ С КОНТРОЛЕМ УСТОЙЧИВОСТИ Если (4.1) выполняется только при l 0, то для построения неравенст ва для контроля устойчивости следует увеличить число вычислений правой части на шаге, т. е. нужно увеличить число стадий m в форму ле (3.19). В этом случае возможно появление свободных коэффициен тов, которые могут быть использованы для улучшения свойств точно сти и устойчивости численной схемы.

Заметим, что проверка неравенства для контроля устойчивости всегда осуществляется после принятия решения по точности, так как в случае невыполнения точности, независимо от устойчивости, величину шага ин тегрирования следует уменьшить. Поэтому применение в (4.1) значения f ( yn+1 ) к увеличению вычислительных затрат не приводит – оно будет использовано при выполнении следующего шага интегрирования.

4.1. СХЕМЫ ВТОРОГО ПОРЯДКА ТОЧНОСТИ Теперь перейдем к получению неравенства для контроля устойчи вости конкретных численных схем. Сначала рассмотрим двухстадий ную формулу yn+1 = yn + (k1 + 3k2 ) / 4, k1 = hf ( yn ), k2 = hf ( yn + 2k1 / 3), (4.2) коэффициенты которой однозначно определяются из требования вто рого порядка точности и условия выполнения второго равенства (2.37), т. е. локальная ошибка численной формулы (4.2) имеет вид n,2 = h3 f 2 f / 6 + O(h3 ). Тогда для контроля точности (4.2) можно применять неравенство || k2 k1 || 4.

Запишем соотношение c1k1 + c2 k2 + c3hf ( yn+1 ) = (c1 + c2 + c3 )hf n + (2c2 / 3 + c3 )h 2 f n f n ) + + h3[c3 f n 2 f n + (4c2 / 9 + c3 ) f n f n2 ] / 2 + O(h 4 ) с неопределенными коэффициентами ci, 1 i 3. Здесь ki, 1 i 2, и yn+1 определены по формуле (4.2). Для выполнения второго условия (4.1) при l = 1 параметры ci, 1 i 3, должны удовлетворять системе уравнений c1 + c2 + c3 = 0, 2c2 / 3 + c3 = 0, 4c2 / 9 + c3 = 0, 4.1. Схемы второго порядка точности которая, очевидно, совместной не является. Таким образом, построить алгоритм интегрирования с контролем точности и устойчивости на ос нове двухстадийной схемы типа Рунге–Кутты нельзя.

Рассмотрим трехстадийную схему вида yn +1 = yn + p1k1 + p2 k2 + p3k3, k1 = hf ( yn ), (4.3) k2 = hf ( yn + 21k1 ), k3 = hf ( yn + 31k1 + 32 k2 ).



Pages:     | 1 || 3 | 4 |   ...   | 9 |
 





 
© 2013 www.libed.ru - «Бесплатная библиотека научно-практических конференций»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.