ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2008, том 48, < 11, с. 2008-2018
УДК 519.622
ЭКОНОМИЧНАЯ СХЕМА РЕАЛИЗАЦИИ НЕЯВНЫХ МЕТОДОВ РУНГЕ-КУТТЫ
© 2008 г. Л. М. Скворцов
(107005 Москва, ул. 2-я Бауманская, 5, МГТУ им. Н.Э. Баумана) e-mail: lm_skvo@rambler.ru Поступила в редакцию 12.04.2008 г.
Предложена схема решения нелинейных алгебраических уравнений при реализации неявных методов Рунге-Кутты. В отличие от известных схем, прогноз стартовых значений для итераций Ньютона выполняется не только для переменных, но и для производных. Это позволяет сократить число вычислений функции (правой части) на каждой неявной стадии до одного без заметного снижения точности интегрирования. Библ. 25. Табл. 11.
Ключевые слова: неявные методы Рунге-Кутты, модифицированные итерации Ньютона, стартовые значения, жесткость.
ВВЕДЕНИЕ
Эффективность решения задачи Коши определяется не только выбранным методом, но и его программной реализацией. В первую очередь это относится к неявным методам, схема реализации которых не задается однозначно коэффициентами формулы интегрирования. Выполнение одного шага неявного метода сводится к численному решению системы нелинейных алгебраических уравнений, для этого обычно применяются итерации Ньютона. Вычислительная схема должна также содержать критерии, согласно которым производится перерасчет матрицы Якоби и прекращение итераций. Вопросы реализации неявных методов Рунге-Кутты рассматривались в работах [1]—[16].
Вычислительные затраты неявного метода можно оценить следующими величинами: 1) число вычислений матрицы Якоби; 2) число вычислений функции; 3) число LU-факторизаций матрицы; 4) число решений линейных уравнений с факторизованной матрицей коэффициентов. Для сложных нелинейных задач наибольшие затраты связаны с вычислением якобиана, поэтому решение алгебраических уравнений производят модифицированным методом Ньютона, в котором матрица Якоби "заморожена" в течение нескольких шагов интегрирования. В этом случае процедура LU-факторизации выполняется не более одного раза на шаге интегрирования (после изменения величины шага или после перерасчета якобиана). Второй по трудоемкости обычно является процедура вычисления функции. Число вычислений функции определяется числом итераций и существенно зависит от выбранного стартового значения. Для сокращения числа итераций применяют специальную формулу прогноза (нетривиальный предиктор, см. [1]—[12]). Для формирования критериев прекращения итераций и обновления матрицы Якоби используют оценку скорости сходимости, полученную на двух последних итерациях (см. [1]).
Таким образом, эффективная схема реализации неявного метода Рунге-Кутты должна использовать замороженную матрицу Якоби и нетривиальный прогноз. Потребуем также, чтобы схема допускала не более одного вычисления функции на каждой неявной стадии, не ухудшая при этом заметно свойств точности и устойчивости исходного метода.
Основная идея предложенной схемы заключается в том, что прогноз выполняется не только для переменных, но и для производных. Формулу прогноза можно получить как значение интерполяционного полинома, построенного по уже вычисленным стадийным значениям (могут быть использованы значения, полученные на предыдущих шагах). В этом случае прогноз для переменных представляет собой линейную комбинацию предыдущих значений переменных. Прогноз для производных получаем в виде такой же линейной комбинации предыдущих значений производных. После окончания итераций определение производных выполняется из формулы интегрирования и не требует вычисления функции.
2008
В предложенной схеме число вычислений функции на каждой неявной стадии на единицу меньше числа итераций. При нетривиальном прогнозе для сходимости бывает достаточно двух итераций, т.е. одного вычисления функции на каждой неявной стадии. Покажем также, что можно построить работоспособную схему, в которой одна из неявных стадий вообще не требует вычислений функции. Рассмотрим подобные схемы на примерах реализации конкретных методов Рунге-Кутты.
1. МЕТОД ТРАПЕЦИИ Будем решать задачу Коши для системы ОДУ
У'= / (г, у), у (го) = Уо, (1.1)
где у - вектор переменных, / - векторная функция, г - независимая переменная. Один шаг решения задачи (1.1) методом трапеций задается формулой
Н
Уп+1 = Уп+2(/п+/п+1 ), (1.2а)
/„ + 1 = /(гп + 1, Уп + 1), (1.2б)
которая представляет собой систему нелинейных алгебраических уравнений относительно Уп + Итерационный процесс вычисления этого вектора запишется в виде
k k -1 f т h Л 1
yn + 1 = yn + 1 + I IJ
2 ( fn + fk + 1 ) - ( yk + 1 - yn )
k = 1, 2,..., N, (1.3a)
kk
fn + 1 = f(tn + 1, Уп + 1), k = 1, 2, ..., N-1, (1.36)
где I - единичная матрица, J - аппроксимация матрицы Якоби, N - число итераций. Будем предполагать, что матрица J сохраняется постоянной в течение нескольких шагов интегрирования и не пересчитывается в процессе итераций.
Определим стартовые значения для итераций. Проще всего задать
yn + 1 = yn, fn + 1 = fn (1.4)
(тривиальный прогноз). Приняв также
N
yk +1 = yk + 1, fn + 1 = f( tn +1, yk +1 ) , (1.5)
получим простейшую схему (1.3)-(1.5), которую назовем тривиальной.
Можно уменьшить число итераций, если использовать нетривиальный прогноз. Применяя экстраполяцию первого порядка, получаем
уП+1 = Уп + W ( Уп + 1- Уп ) , f П + 1 = f (tn + 1, У П + 1 ) , (1.6)
где w = h/h - отношение размера текущего шага к размеру предыдущего шага. Вектор производных fk + 1 будем определять из формулы интегрирования (1.2а), а не из (1.26), тогда
2
Уп +1 = yf+1, fn + 1 = J- ( Уп + 1- Уп ) - fn • (1.7)
Такой прием применяется при реализации неявных методов с явной первой стадией; в [13] он назван "сглаженная первая стадия" (smoothed first stage). Схему, задаваемую формулами (1.3), (1.6), (1.7), назовем стандартной. Подобные схемы обычно используются в программных реализациях неявных методов Рунге-Кутты. В стандартной схеме, как и в тривиальной, число вычислений функции равно числу итераций.
В предлагаемой схеме прогноз выполняется не только для переменных, но и для производных, т.е. вместо (1.6) принимаем
УП+ 1 = Уп + w(yn +1-Уп ), fk+1 = fn + W( fn + 1- fn) • (1.8)
Схема
10° 101 102 103 104 105
T(1) 7.0 x 10-3 2.2 x 10-2 2.8 x 10-2 2.9 x 10-2 2.9 x 10-2 2.9 x 10-2
T(2) 2.1 x 10-4 2.8 x 10-4 2.1 x 10-4 7.9 x 10-5 5.0 x 10-5 4.7 x 10-5
T(3) 1.5 x 10-4 7.5 x 10-5 5.0 x 10-5 4.7 x 10-5 4.6 x 10-5 4.6 x 10-5
S(1) 1.7 x 10-4 6.8 x 10-4 9.0 x 10-4 9.2 x 10-4 9.1 x 10-4 9.2 x 10-4
S(2) 1.5 x 10-4 6.6 x 10-5 4.4 x 10-5 4.6 x 10-5 4.6 x 10-5 4.6 x 10-5
E(1) 1.5 x 10-4 7.0 x 10-5 4.7 x 10-5 4.6 x 10-5 4.6 x 10-5 4.6 x 10-5
Точная 1.5 x 10-4 7.3 x 10-5 4.9 x 10-5 4.7 x 10-5 4.6 x 10-5 4.6 x 10-5
Схему, задаваемую формулами (1.3), (1.8), (1.7), назовем экономичной, в ней число вычислений функции на единицу меньше числа итераций.
Для тестирования рассмотренных схем использовалась задача Капса
У1 = -(Ц + 2) У1 + Ц У22, 0) = 1,
У2 = У1-У2-У22, У2(0) = 1, 0 < г < 1,
имеющая гладкое решение У1(г) = ехр(-2г), У2(г) = ехр(-г), не зависящее от параметра жесткости ц. На этой задаче удобно исследовать зависимость ошибки от жесткости. Определялась максимальная относительная ошибка на всем интервале интегрирования при условии, что матрица Якоби вычисляется только один раз в начальной точке. Ошибки решения задачи Капса различными схемами метода трапеций с шагом Н = 1/60 приведены в табл. 1, где приняты обозначения: Т - тривиальная схема, Б - стандартная схема, Е - экономичная схема (в скобках указано число вычислений функции на одном шаге интегрирования). В последней строке приведены ошибки при точной реализации метода трапеций, т.е. при выполнении итераций до сходимости. Такие ошибки получены уже при использовании схем Т(4), Б(3), Е(2). В этой и последующих таблицах максимальная ошибка в каждой строке выделена жирным шрифтом.
Из приведенных результатов видно, что сходимость итераций экономичной схемы практически обеспечивается при одном вычислении функции, тогда как тривиальная схема требует трех вычислений, а стандартная схема - двух вычислений. Отметим, что самая медленная сходимость наблюдается при умеренной жесткости (ц = 101, 102).
2. ДИАГОНАЛЬНО НЕЯВНЫЕ МЕТОДЫ
Наиболее просто реализуются однократно диагонально неявные методы Рунге-Кутты, среди которых преимущество имеют жестко точные методы ЕЗБЖК второго стадийного порядка (см. [13]—[18]). Формулы одного шага метода ЕЗБШК можно представить в виде
У 0 = Уп, = /п, (2.1а)
Yi = Уп + h
X aijFj + Y Fi
Ч = 0
Fi = f (tn + chi, Yi), i = 1, 2, ..., 5,
yn + 1 Y5, fn + 1 F 5,
(2.16)
(2.1B)
где Yi и Fi - векторы стадийных значений переменных и производных. Такой метод имеет одну явную стадию (2.1а) и s неявных стадий (2.16). Он является жестко точным, что выражается заключительной формулой (2.1в). Результат выполнения явной стадии (2.1а) совпадает с результатом выполнения последней стадии предыдущего шага, поэтому методы такого типа называют также FSAL (First stage Same As Last). Для удобства изложения присваиваем явной стадии индекс 0
(нулевая стадия). Обозначим через Yi и Fi стадийные значения на предыдущем шаге.
ЭКОНОМИЧНАЯ СХЕМА РЕАЛИЗАЦИИ НЕЯВНЫХ МЕТОДОВ РУНГЕ-КУТТЫ 2011 Целесообразно вместо стадийных значений использовать приращения
г{ = У1 - Уп, I = 1, 2, ..., 5; (2.2) тогда итерации Ньютона на г-й стадии запишутся в виде
гк = гк-1 + (I - ум)-
-1
Н
X р + У Г
\р = о
л-1
- г
к-1
, к = 1, 2,..., и,
7к = Уп + гк, Гк = /(гп + с,Н, 7к), к = 1, 2, ..., N-1. В стандартной схеме в качестве стартовых значений принимаем
5 -1 г -1
7° = + ХМ,
р = °
р = °
Г° = /(гп + ф, 7°),
(2.3а) (2.36)
(2.4а) (2.46)
где некоторые коэффициенты ар, вр могут быть нулевыми. Завершается стадия вычислением стадийных значений
N 1
г = г- г = 1-
г-1
Н г - X р
р = 0 -
(2.5)
Формулы (2.3)-(2.5) задают стандартную схему вычислений на одной неявной стадии. В экономичной схеме вместо (2.46) будем использовать
5-1 г-1
Г° = Хр+ХврГ.
р=0 р=0
Пусть 7* - точное решение, к которому сходятся итерации (2.3). Согласно [3], [4], [9], прогноз
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.