CC BY
44
1. Харебов, К.С. Компьютерные методы решения задачи наименьших квадратов и проблемы собственных значений / К.С. Харебов. - Владикавказ: Изд-во СОГУ, 1995. - 76 с.
2. Лоусон, Ч. Численное решение задач наименьших квадратов / Ч. Лоусон, Р. Хенсон. - М.: Статистика, 1989. - 447 с.
3. Буераков, Н.Я. Организация финансирования лесного хозяйства за счет средств лесного дохода (на примере Удмуртской Республики): дис. ... канд. экон. наук / Н.Я. Буераков. - М., 1995.
4. Бурдин, Н.А. Лесопромышленный комплекс. Состояние, проблемы, перспективы / Н.А. Бурдин, В.М. Шлыков, В.А. Егорнов, В.В. Саханов. - М.: МГУЛ, 2000. - 473 с.
5. Воронин, И.В. Экономика лесного хозяйства / И.В. Воронин, П.В. Васильев, Е.Я. Судачков. - М.: Лесная промышленность, 1978. - 264 с.
6. Ильин, В.А. Вопросы экономической организации промежуточного пользования лесом: автореф. дис. ... канд. экон. наук / В.А. Ильин. - Л., 1973.
УДК 519.622 Е.А. Новиков, Л.В. Кнауб
АЛГОРИТМ ИНТЕГРИРОВАНИЯ НА ОСНОВЕ ЯВНОГО ТРЕХСТАДИЙНОГО МЕТОДА РУНГЕ-КУТТА*
Разработан алгоритм интегрирования переменного шага на основе трехстадийной схемы типа Рунге-Кутта. Численная схема имеет третий порядок точности. Получена оценка аналога глобальной ошибки, которая не приводит к увеличению вычислительных затрат. На основе данной оценки построено неравенство для контроля точности вычислений и автоматического выбора величины шага интегрирования. Сформулирован алгоритм управления величиной шага интегрирования исходя из требования точности.
Ключевые слова: методы Рунге-Кутта, оценка глобальной ошибки, контроль точности вычислений.
E.A. Novikov, L.V. Knaub INTEGRATING ALGORITHM BASED ON THE EXPLICIT THREE- STAGE RUNGE-KUTTA METHOD
Integrating algorithm of a variable step integration based on the explicit three- stage Runge-Kutta scheme is developed. The numerical scheme has the third order of accuracy. The estimation of the analogy which does not increase calculating expenses is obtained. On the basis of the estimation the inequality for the accuracy control of the calculations and for the automatic choice of the integrating step size is developed. The algorithm for the integrating step size control is formulated on the basis of the accuracy requirement.
Key words: Runge-Kutta methods, global error appraisal, calculations exactness control.
При решении задачи Коши для жестких систем обыкновенных дифференциальных уравнений большой размерности в ряде случаев возникает необходимость применения алгоритмов на основе явных методов [1-3]. Алгоритмы интегрирования на основе неявных или полуявных формул, как правило, используют обращение матрицы Якоби [2], что в данном случае есть отдельная трудоемкая задача. Кроме того, получение элементов матрицы Якоби и составление подпрограммы ее нахождения требуют от вычислителя больших затрат личного времени. В такой ситуации предпочтительнее применять алгоритмы на основе явных формул, если жесткость задачи позволяет за разумное время получить приближение к решению [3]. Обычно алгоритм управления шагом интегрирования строится на контроле точности численной схемы. Это естественно, потому что основным критерием является точность нахождения решения.
В данной работе разработан алгоритм интегрирования переменного шага на основе трехстадийной схемы типа Рунге-Кутта. Численная схема имеет третий порядок точности. Получена оценка аналога глобальной ошибки, которая не приводит к увеличению вычислительных затрат. На основе данной оценки построено неравенство для контроля точности вычислений и автоматического выбора величины шага интегрирования. Сформулирован алгоритм управления величиной шага интегрирования исходя из требования точности.
* Работа поддержана грантами РФФИ №08-01-00621 и Президента НШ-3431.2008.9.
Для численного решения задачи Коши для системы обыкновенных дифференциальных уравнений
/= , = , < < (D
рассмотрим явный трехстадийный метод типа Рунге-Кутта вида
Уп+ = .+ . . + . . + _
К= _= ..+ .» (^)
к^ = + + ,
где у и f - вещественные N-мерные вектор-функции;
t - независимая переменная; h - шаг интегрирования; кх, к2 и къ ~ стадии метода;
рх, , , , и (3 - числовые коэффициенты, определяющие свойства точности и устойчиво-
сти (2).
В случае неавтономной системы
У = , ,
схема (2) записывается в виде
< <
Уп+ - + . . + . . + . .
К = _= ..+ > ..+
кз — + + , + + .
Ниже для сокращения выкладок будем рассматривать задачу (1). Однако построенные методы можно применять для решения неавтономных задач.
Получим соотношения на коэффициенты метода (2) третьего порядка точности. Для этого разложим
стадии £, к2 и £ в ряды Тейлора до членов с к4 включительно, т.е.
К = >
л л
— + + + + ,
" 2 “* " ■' Ь - " "
— + + + +
1
+ + " . +
2 " '■ " ” _ '■ ’ " " "
1
+ + + + +0(И5).
' -- - - ......... 5 - ......
Подставим полученные разложения в первую формулу (2). В результате получим
Уп+ = . + . + _ + _ .. + .. .+ _ .. .. +
+ + . + + + (3)
+ . + + +
где элементарные дифференциалы вычислены на приближенном решении уп, /п=д /д , /„=д /д , /и = д /д ■ Точное решение у(ги+ 4 в окрестности точки (п имеет вид
1
7.3
ж+ч= .. + + • + ^ . + +
ь
1А
+ +++ + , (4)
24
где элементарные дифференциалы вычислены на точном решении у(гп). Сравнивая полученные ряды для точного и приближенного решений до членов с /г3 включительно при условии уп = , запи-
шем условия третьего порядка точности схемы (2), которые имеют вид
1
Л+ _+ -= > __ _ + -=, 5
р + + =1, = 1. (5)
■ ■ “ 3 ()
Теперь локальную ошибку § схемы (2) можно вычислить по формуле
8 = - .
Учитывая представления у(ґп+ Л и уп+ в виде рядов Тейлора, получим
1 11
5 = + - +
" 24' 24 2 “ ' '
1
+ - + + (6)
8 - -- - - -
її і + - - + + .
24 6 “ “ Ь "
В нелинейной системе алгебраических уравнений (2) имеются два свободных коэффициента. Исследуем три варианта.
Вариант 1. Положим р = + и р = . Это означает, что приращения к2 и къ будут вычислены в одной и той же точке ^ + , причем вклад кх и к2 при определении к, учитывается одинако-
во. Тогда нелинейную систему (5) можно переписать в виде
1
1)# + _ + _ = >
3)£_ _ + 1
з о
Из второго и третьего уравнений имеем р = / . Из соотношений р = + и р = запишем р = = / . Из четвертого уравнения системы получим р3= / . Из равенства
р2+ = / имеем р2 = . Наконец, из первого соотношения системы получим рх = / .В результате
коэффициенты схемы (2) определяются однозначно и имеют вид
п 11 о
р = , = = , = , = , = . (7)
- 3 " - '3 ‘4 ■ '4
При данных соотношениях локальную ошибку 8 схемы (2) можно записать следующим образом
Л Л л л
5 = + - + + ■
24 210 72 72
Вариант 2. Минимизируем локальную ошибку (6). Для этого, учитывая вид (6), вместо (5) рассмотрим следующую расширенную нелинейную систему:
1
1)А + _+ _= > _ + .. + __ ~= '2>
3)/? + + = , = , (8)
... .. .. . 3 .. .. . ь
5)0 (3 = 1 , + = ’ .
__ . 12 __ .. „ „ . 8
При 1,5[3 , = 3 , + 3 2 два последних уравнения (8) совпадают. Из четвертого и пятого соотношений (8) имеем Р = ),5. Из второго и третьего равенств получим р2= / и р3= / .Из первого уравнения (8) запишем р] = / , а из четвертого имеем р = / . Наконец, из соотношения р + = / запишем р = .
В результате коэффициенты метода (2) с минимальной локальной ошибкой можно записать в виде
1 0 0 1/1 Р = , = , =,=,=,=. (9)
2 4 У "3 У
При данных соотношениях локальную ошибку 8 схемы (2) можно представить следующим образом:
7.4 7 4
3 = - + ■
24 28»
При использовании (2) с наборами коэффициентов (7) и (9) ни одна стадия не вычисляется в точке tn+ . При быстром изменении решения это может приводить к понижению точности расчетов.
Вариант 3. Положим р = . и р + = . Тогда на каждом шаге приращения кх, к2 и къ вычисля-
ются в точках tn,t + . я (п+ , соответственно. В этом случае условия третьего порядка записываются в виде
1 1
1 )Р\ + + = 5 + = 5
1 . . 2 ” ' 2
1 1 1
3 )-лР2+ = , = . (10)
4 ' 3 - - 3
Из второго и третьего уравнения данной системы имеем р = / и ръ = / . Из первого и последнего уравнений имеем рх = ! и р = .Из равенства /3 + = следует р = - . В результате
коэффициенты метода (2) можно записать в следующем виде:
р = =- = =’ =" =1 (11)
? Г> 9 9 ? *
“ 2 " Ь ~ 5 Ь
При данных соотношениях локальную ошибку 8 схемы (2) можно записать следующим образом:
7.4 ,
5 = - - + ■
24 3
Построим неравенство для контроля точности вычислений метода третьего порядка. Для этого рассмотрим вспомогательную схему
Уп+. =
(12)
где кх и к2 определены в (2). Потребуем, чтобы данный метод имел второй порядок точности. Разложение
приближенного решения уп+ в виде ряда Тейлора по степеням к имеет вид
Уп+ , =
Сравнивая ряды для у(/п+ ' и уп видим, что данное требование будет выполнено, если
Г1 + _ = ’ _ = • •
Отсюда получим
г2= . / (13)
где значение р определено в (7), (9) или (11).
С помощью идеи вложенных методов оценку ошибки £ метода третьего порядка можно оценить по формуле
<14)
Тогда неравенство для контроля точности вычислений имеет вид
II (А - . . + ^ , (15)
где |||| - некоторая норма в , £ - требуемая точность интегрирования. В конкретных расчетах применялся метод (2) с коэффициентами (11), как наиболее надежный. Тогда неравенство для контроля точности имеет вид
-II*!- + ^ , (16)
6 '
Сформулируем алгоритм интегрирования переменного шага на основе схемы (2). Заметим, что
£ =0(№), данный факт будет использоваться при выборе шага интегрирования.
Шаг 1. Вычисляется стадия Ь по формуле (2).
Шаг 2. Вычисляются стадии к2 и кз по формулам (2).
Шаг 3. Вычисляется оценка ошибки еп,з по формуле (14).
Шаг 4. Вычисляется число q по формуле
с£е = ,
где г - требуемая точность интегрирования.
Шаг 5. Если q<1, то есть требуемая точность не выполнена, то Ь полагается равным qh и происходит возврат на шаг 2.
Шаг 6. Вычисляется приближенное решение по первой формуле (2), величина шага Ь полагается равной qh.
Шаг 7. Выполняется следующий шаг интегрирования - переход на шаг 1.
Литература
1. Хайрер, Э. Решение обыкновенных дифференциальных уравнений. Нежесткие задачи / Э. Хайрер, С. Нерсетт, Г. Ваннер. - М.: Мир, 1990. - 512 с.
2. Хайрер, Э. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциальноалгебраические задачи / Э. Хайрер, Г. Ваннер. - М.: Мир, 1999. - 685 с.
3. Новиков, Е.А. Явные методы для жестких систем / Е.А. Новиков. - Новосибирск: Наука, 1997. - 197 с.
