close

Вход

Забыли?

вход по аккаунту

?

Интегратор Гаусса-Эверхарта.

код для вставкиСкачать
Вычислительные технологии
Том 15, ќ 4, 2010
Интегратор ауссаЭверхарта
?
В. А. Авдюшев
НИИ прикладной математики и механики
Томского государственного университета, оссия
e-mail: shniipmm.tsu.ru
Представлена новая версия интегратора Эверхарта для решения обыкновенных
диеренциальных уравнений первого порядка. Обсуждаются особенности в компьютерной реализации интегратора, исследуются его возможности при решении
задач небесной механики.
Ключевые слова: численное интегрирование, обыкновенные диеренциальные
уравнения, интегратор Эверхарта, неявные методы унгеКутты, численное моделирование орбит.
Введение
В 1973 г. Э. Эверхарт [1? предложил интегратор, специально разработанный для численного исследования орбит, и продемонстрировал его высокую эективность в задачах
кометной динамики. По-видимому, обнаружив в дальнейшем принадлежность данного интегратора к семейству интеграторов бутчеровского типа, Эверхарт акцентировал
внимание на оригинально реализованный им алгоритм интегрирования и обобщил его
для численного решения любых обыкновенных диеренциальных уравнений первого и второго порядка [2, 3?, тем самым расширив область применения предложенного
интегратора, который, тем не менее, остается одним из самых популярных именно в
решении задач небесной механики.
Интегратор Эверхарта (RA15) основан на видоизмененных ормулах неявных коллокационных методов унгеКутты бутчеровского типа, поэтому он наследует все их
замечательные свойства [4?. Более того, именно благодаря оригинальному представлению вычислительной схемы интегратор Эверхарта с точки зрения численного интегрирования имеет ряд следующих преимуществ:
1) алгоритм интегрирования универсален для любого порядка;
2) интегратор имеет простой критерий для выбора шага интегрирования;
3) в интеграторе реализован достаточно точный предиктор решения, что позволяет
выполнять численное интегрирование всего с двумя итерациями на шаге.
Несмотря на это, программный код Эверхарта RA15 [3? (впрочем, как и любая
его модиикация типа RADAU_27), на наш взгляд, довольно существенно ограничивает возможности интегратора и поэтому нуждается в дополнительной редакции.
?
абота выполнена при инансовой поддержке ФФИ (грант ќ 08-02-00359) и Министерства об-
разования Ф (грант ќ НП.2.1.1/2629).
31
32
В. А. Авдюшев
Среди главных недостатков в программной реализации интегратора можно назвать
следующие:
1) трудночитаемый и громоздкий код;
2) много констант, связанных с порядком интегратора, что затрудняет обобщение
кода на другие порядки;
3) интегратор реализован только для определенных порядков, причем для нечетных с разбиением ауссаадо, хотя известно, что неявные коллокационные методы
унгеКутты, построенные на симметричных разбиениях ауссаЛобатто и аусса
Лежандра, обладают геометрическими свойствами [5?;
4) алгоритм выбора шага для уравнений первого порядка используется такой же,
как и для уравнений второго порядка, поэтому шаг при интегрировании уравнений
первого порядка выбирается неверно;
5) стартовый шаг интегрирования в режиме переменного шага выбирается независимо от диеренциальных уравнений, поэтому не всегда оптимально;
6) ограничения на величину выбираемого переменного шага, по-видимому, заданы в
интеграторе просто из эмпирических соображений и, кроме того, не зависят от порядка
интегратора.
В настоящей работе представлен новый код интегратора Эверхарта GAUSS_15, который позволяет устранить перечисленные выше недостатки:
1) путем использования возможностей Фортран 90 программный код сокращен почти в 2 раза;
2) устранены все константы, связанные с порядком метода (оставлены лишь константы узловых значений на шаге);
3) код позволяет получать решение 215 порядка точности (хотя при необходимости
код без изменений можно обобщить на любой другой порядок: для этого нужно лишь
получить соответствующие узловые значения);
4) исправлен алгоритм выбора переменного шага;
5) стартовый шаг выбирается по оценке интегрирующей схемы второго порядка с
учетом поведения правых частей уравнений;
6) накладываются ограничения на выбираемый шаг в соответствии с порядком интегратора.
Кроме того, предложенный интегратор имеет новые возможности:
1) интегрирование на шаге до полной сходимости итерационного процесса;
2) запоминание величины предпоследнего шага после выполнения процедуры интегрирования при многократном использовании программного кода в режиме переменного шага;
3) быстрый выбор стартового шага, требуемый лишь для первого обращения к интегратору (при повторном обращении применяется запоминаемый шаг предыдущего
обращения).
В работе представлена общая теория интегратора Эверхарта с внесенными автором
коррективами, предлагается новый программный код интегратора, а также показаны
результаты тестирования интегратора на примере диеренциальных уравнений задачи двух тел. В дальнейшем ввиду того что интегратор использует гауссовы разбиения,
будем называть его интегратором ауссаЭверхарта (хотя обычно такие интеграторы
называются гауссовыми [5?).
Интегратор ауссаЭверхарта
33
1. Основные ормулы
Предположим, что на шаге h мы решаем задачу
x? = f(t, x),
x0 = x(t0 ).
(1)
Здесь t независимая переменная, x интегрируемые переменные, f заданная
векторункция t и x. Введем переменную ? = (t ? t0 )/h и представим правую часть
уравнений (1) в виде полинома степени k :
x? = x?? /h = f = f0 + A1 ? + A2 ? 2 + A3 ? 3 + . . . + Ak ? k ,
(2)
где коэициенты Ai пока не определены. Интегрируя (2) по ? , получаем решение
1
1
1
1
2
3
4
k+1
x = x0 + h f0 ? + A1 ? + A2 ? + A3 ? + . . . +
Ak ?
.
(3)
2
3
4
k+1
Перепишем (2) в виде интерполяционного многочлена Ньютона на сетке ?0 , ?1 , . . . , ?k
(?0 = 0):
f = f0 + ?1 ? + ?2 ? (? ? ?1 ) + ?3 ? (? ? ?1 )(? ? ?2 ) + . . . + ?k ? (? ? ?1 ) . . . (? ? ?k?1 ). (4)
Из соотношений
f1 = f0 + ?1 ?1 ,
f2 = f0 + ?1 ?2 + ?2 ?2 (?2 ? ?1 ),
f3 = f0 + ?1 ?3 + ?2 ?3 (?3 ? ?1 ) + ?3 ?3 (?3 ? ?1 )(?3 ? ?2 ),
...
(5)
получаем конечные разности ?
?1 = (f1 ? f0 )/?1 ,
?2 = ((f2 ? f0 )/?2 ? ?1 )/(?2 ? ?1 ),
?3 = (((f3 ? f0 )/?3 ? ?1 )/(?3 ? ?1 ) ? ?2 )/(?3 ? ?2 ),
...
(6)
В свою очередь, сравнивая (2) и (4), будем иметь
A1 = ?1 + (??1 )?2 + (?1 ?2 )?3 + . . . + (?1)k?1 (?1 . . . ?k?1 )?k ,
A2 = ?2 + (??1 ? ?2 )?3 + . . . ,
...
A k = ?k ,
или
A1 = c11 ?1 + c21 ?2 + . . . + ck1 ?k ,
A2 = c22 ?2 + . . . + ck2 ?k ,
...
Ak = ckk ?k .
(7)
В. А. Авдюшев
34
Тогда обратный переход от A к ? можно представить как
?1 = d11 A1 + d21 A2 + . . . + dk1 Ak ,
?2 = d22 A2 + . . . + dk2Ak ,
...
?k = dkk Ak .
(8)
Коэициенты cij и dij являются числами Стирлинга и вычисляются по ормулам
cii = dii = 1, ci0 = di0 = 0 (i > 0),
cij = ci?1,j?1 ? ?i?1 ci?1,j , dij = di?1,j?1 ? ?j di?1,j
(i > j > 0).
(9)
В соответствии с (9) каждый коэициент cij представляет собой сумму всевозможных
произведений i ? j величин ?1 , . . . , ?i?1 со знаком (?1)i?j , например,
c62 = ?1 ?2 ?3 ?4 + ?1 ?2 ?3 ?5 + ?1 ?2 ?4 ?5 + ?1 ?3 ?4 ?5 + ?2 ?3 ?4 ?5 .
Отсюда согласно теореме Виета значения ?1 , . . . , ?i?1 будут являться корнями полинома
вида
Pi?1 = ci1 + ci2 ? + ci3 ? 2 + . . . + cii ? i?1 .
(10)
2. Интегрирование на шаге
Величины ? определяются по f , которые, в свою очередь, вычисляются по решениям
x. Согласно (3) эти решения будем задавать как
1
1
1
2
3
k+1
Ak ?1
,
x1 = x0 + h f0 ?1 + A1 ?1 + A2 ?1 + . . . +
2
3
k+1
...
1
1
1
2
3
k+1
xk = x0 + h f0 ?k + A1 ?k + A2 ?k + . . . +
Ak ?k
.
(11)
2
3
k+1
Формулы (11) представляют собой неявные уравнения относительно x, поэтому решаются итерационным способом.
Для получения начального приближения ?? на следующем шаге h? используется инормация о коэициентах A на текущем шаге h. Безразмерная независимая переменная следующего шага будет ?? = (t ? th )/h?, где th = t0 + h. Отсюда
? = r?? + 1,
(12)
где r = h?/h. Согласно (2)
f0 + A1 ? + A2 ? 2 + A3 ? 3 + . . . + Ak ? k = f?0 + A?1 ?? + A?2 ?? 2 + A?3 ?? 3 + . . . + A?k ?? k .
(13)
Подставляя (12) в (13) и приравнивая коэициенты при одинаковых степенях ?? , получаем
f?0 = e00 f0 + e10 A1 + e20 A2 + e30 A3 + . . . + ek0 Ak ,
A?1 = r(e11 A1 + e21 A2 + e31 A3 + . . . + ek1 Ak ),
A?2 = r 2 (e22 A2 + e32 A3 + . . . + ek2 Ak ),
...
A?k = r k ekk Ak ,
(14)
Интегратор ауссаЭверхарта
35
где eij числа ариметического треугольника, вычисляемые по рекуррентным ормулам
eii = ei0 = 1, eij = ei?1,j?1 + ei?1,j (i > j > 0).
(15)
Далее оценка A? для h? уточняется путем внесения поправки ?A, получаемой как разность между значениями A после итераций и оценкой A? на текущем шаге h. Наконец,
пользуясь соотношениями (8), будем иметь начальное приближение ??.
Каждая итерация выполняется следующим образом. Сначала определяется решение x1 , из которого по первой ормуле (5) улучшается значение ?1 . Далее определяется x2 , по которому улучшается ?2 , и т. д. до xk . Как правило, для получения достаточно
хороших ? необходимы всего две итерации, очень редко три.
Как только величины ? получены, решение на шаге h (? = 1) будет
1
1
1
xh = x0 + h f0 + A1 + A2 + . . . +
Ak .
2
3
k+1
(16)
В начале интегрирования, на первом шаге, в качестве ?? выбирают нулевые значения
и запускается вышеописанный итерационный процесс. Если начальный шаг достаточно
большой, чтобы обеспечить заданную локальную точность, то его следует уменьшить.
При оптимально выбранном шаге высокая точность ? достигается уже на четвертой
итерации.
3. Формулы интегратора как одно из представлений
неявного метода унгеКутты
Пользуясь (6), (7) и (11), нетрудно показать, что решение (16) представимо в виде
xh = x0 + h
k
X
bi ki ,
где
ki = f(t0 + h?i , x0 + h
i=0
k
X
aij kj ) (i = 0, . . . , k),
(17)
j=0
а коэициенты aij и bi постоянные, зависящие только от ?i . Таким образом, интегратор ауссаЭверхарта актически основан на видоизмененных ормулах неявного
метода унгеКутты, причем эти ормулы являются коллокационными [6?. Действительно, в качестве коллокационного полинома в интеграторе выступает приближенное
представление решения (3), тогда как условия коллокации ормально задаются соотношениями (5).
Произвольный выбор узлов ?i дает метод порядка k + 1. Однако, используя свойства
неявных методов унгеКутты, Эверхарт выбрал разбиение ауссаадо, что позволило повысить точность метода до порядка 2k + 1.
4. Повышение порядка интегратора
Значения узлов ?1 , . . . , ?k свободные параметры и их можно выбирать как угодно,
лишь бы они не равнялись между собой. В общем случае интегратор будет иметь порядок k + 1, однако существуют такие узловые значения, которые позволяют значительно
повысить точность решения до порядка 2k либо 2k + 1.
В. А. Авдюшев
36
Приближенное решение порядка 2k + 1 на сетке ?1 , . . . , ?2k будет
1 ?
1
1 ?
?
xh = x0 + h f0 + A1 + A2 + . . . +
A
.
2
3
2k + 1 2k
(18)
Определим ?1 , . . . , ?k так, чтобы решение (16) имело порядок 2k + 1. Это возможно в
случае, если разность решений (16) и (18) будет равна нулю:
?
A?k ? Ak A?k+1
A?2k
A1 ? A1 A?2 ? A2
?xh = h
+
+ ...+
+
+ ...+
= 0.
(19)
2
3
k+1
k+2
2k + 1
Согласно (6) коэициенты ?1 , . . . , ?k в (16) и (18) совпадают. Тогда
A?1 ? A1 = ck+1,1?k+1 + . . . + c2k,1 ?2k ,
A?2 ? A2 = ck+1,2?k+1 + . . . + c2k,2 ?2k ,
...
A?k ? Ak = ck+1,k ?k+1 + . . . + c2k,k ?2k ,
A?k+1 = ck+1,k+1?k+1 + . . . + c2k,k+1?2k ,
A?k+2 = ck+2,k+2?k+2 + . . . + c2k,k+2?2k ,
...
A?2k = c2k,2k ?2k .
(20)
Используя рекуррентные соотношения (9), выразим в (20) коэициенты ck+i,j через
ck+1,j :
ck+2,1 = ??k+1 ck+1,1 , ck+2,j = ck+1,j?1 ? ?k+1 ck+1,j (j > 1),
ck+3,1 = ?k+2 ?k+1 ck+1,1 , ck+3,j = ck+1,j?2 + (??k+1 ? ?k+2 )ck+1,j?1 + ?k+2 ?k+1 ck+1,j (j > 1),
...
(21)
Подставляя (21) в (20), получим
A?1 ? A1 = B1 ck+1,1 ,
A?2 ? A2 = B1 ck+1,2 + B2 ck+1,1 ,
...
A?k ? Ak = B1 ck+1,k + B2 ck+1,k?1 + . . . + Bk ck+1,1 ,
A?k+1 = B1 ck+1,k+1 + B2 ck+1,k + . . . + Bk ck+1,2 ,
A?k+2 = B2 ck+1,k+1 + . . . + Bk ck+1,3,
...
A?2k = Bk ck+1,k+1,
где B имеют орму (7):
B1 = ?k+1 + (??k+1 )?k+2 + (?k+1 ?k+2 )?k+3 + . . . + (?1)k?1 (?k+1 . . . ?2k?1 )?2k ,
B2 = ?k+2 + (??k+1 ? ?k+2 )?k+3 + . . . ,
...
Bk = ?2k .
(22)
Интегратор ауссаЭверхарта
37
Подставляя далее (22) в (19) и уравнивая коэициенты при одинаковых величинах
B, получим следующие уравнения для ck+1,j :
ck+1,1 ck+1,2
+
+ ...+
2
3
...
ck+1,1 ck+1,2
+
+ ...+
k+1 k+2
ck+1,k
1
+
= 0,
k+1 k+2
ck+1,k
1
+
= 0.
2k
2k + 1
(23)
ешения (23) однозначно определяют значения ?1 , . . . , ?k , которые вычисляются из уравнения
ck+1,1 + ck+1,2 ? + ck+1,3 ? 2 + . . . + ck+1,k ? k?1 + ? k = 0.
(24)
Узловые значения ?i , определяемые из (23) и (24), задают (левое) разбиение аусса
адо. Их также можно получить из уравнения
(25)
(? k+1 (? ? 1)k )(k)
? = 0.
В (25) левый корень равен нулю (?0 = 0). Если потребуем, чтобы правый корень был
равен единице (?k = 1), то для повышения порядка численного решения получим разбиение ауссаЛобатто, узловые значения которого также удовлетворяют уравнению
(26)
(? k (? ? 1)k )(k?1)
= 0.
?
азбиение ауссаЛобатто симметрично и дает численное решение порядка 2k .
5. Выбор шага
В интеграторе ауссаЭверхарта контроль шага интегрирования осуществляется по
величине последнего члена в (16).
Пусть ketol k заданная точность. Потребуем, чтобы на следующем шаге выполнялось равенство
h?
kA?k k = ketol k.
k+1
Отсюда, используя последнее соотношение в (14), получим оценку
h? = hr = h
k + 1 ketol k
h kAk k
1
k+1
.
(27)
Очевидно, при разбиениях ауссаадо и ауссаЛобатто недостаток такой оценки
состоит в том, что шаг по ней выбирается как для решения порядка k , поэтому, вообще
говоря, она не обеспечивает сохранение требуемой локальной точности для решения
более высокого порядка.
Во избежание слишком больших (и малых) локальных ошибок на r следует наложить ограничение:
1
< r k+1 < ?.
(28)
?
Чтобы величина последнего члена
одного порядка,
? в (16) была ограничена в пределах
k+1
значение ? должно быть равно 10. Это следует из h?kA?k k ? r .
В. А. Авдюшев
38
Выполнение обоих неравенств проверяется лишь в начале интегрирования при выборе стартового шага: если (28) не выполняется, то интегрирование повторяется с новым
шагом h? = hr и т. д., пока не выполнится условие (28). Обычно для получения стартового шага требуется не более четырех итераций. В дальнейшем для ограничения r
проверяется только правое неравенство: если оно не выполняется, то r принимает значение правого предела.
Начальное приближение стартового шага получается из оценки (27) для k = 1:
s
2hketol k
h? =
, f1 = f(t0 + h, x0 + hf0 ),
(29)
kf1 ? f0 k
где h малая величина. Если h настолько мала, что в компьютерной ариметике
f1 = f0 , то она увеличивается в 10 раз и оценка повторяется снова.
6. Порядок и шаг интегрирования
при компьютерной реализации метода
Теоретически совместное повышение порядка и уменьшение шага метода неограниченно повышает методическую точность численных результатов интегрирования. Однако при компьютерной реализации в ариметике с определенной точностью вследствие
ошибок округления существуют такие значения параметров интегрирования, которые
дают предельно высокую точность ввиду того, что методические ошибки становятся
соизмеримыми с ошибками округления, и в этом случае не имеет смысла предпринимать какие-либо дальнейшие попытки получить более высокую точность численного
интегрирования путем повышения порядка и уменьшения шага метода.
В общем случае получить оптимальные параметры интегрирования невозможно,
поскольку они непосредственно зависят от специики решаемой задачи. Впрочем, принимая во внимание, что интегратор ауссаЭверхарта используется, как правило, для
численного моделирования орбит, найдем оценки оптимальной пары порядокшаг интегрирования применительно к решению какой-нибудь простой задачи небесной механики, например, круговой двумерной задачи двух тел с уравнениями
r? = v,
v? = ?
µ
r.
|r|3
(30)
Здесь r = (r1 , r2 ), v = (v1 , v2 ) векторы положения и скорости соответственно, µ гравитационный параметр. Поскольку в круговом случае |r| = a = onst, будем полагать,
что величина µ/|r|3 = n2 = onst, т. е. решение x = (r, v) описывается уравнениями
гармонического осциллятора с частотой n и может быть представлено в виде
x1 = r1 = a cos nt,
x3 = v1 = ?an sin nt,
x2 = r2 = a sin nt,
x4 = v2 = an cos nt.
Оценим методическую ошибку |e|M по главному члену погрешности:
?
a 1 + n2 (nh)p+1
|e|M =
,
(p + 1)!
(31)
(32)
Интегратор ауссаЭверхарта
39
?
где использована ормула |x(p+1) | = anp+1 1 + n2 . Здесь p порядок метода. Ошибку
округления |e|R можно оценить как
?
|e|R = ?|x| = ?a 1 + n2 ,
(33)
где ? машинный эпсилон.
Очевидно, что не имеет смысла выбирать такие шаг и порядок интегрирования, при
которых методическая ошибка будет меньше ошибки округления. Из условия |e|M =
|e|R получим отношение между оптимальными параметрами интегрирования h и p:
(34)
hp+1
= ?(p + 1)!,
?
где h? = nh шаг по долготе ? = nt, соответствующий шагу h. Уравнение (34) дает
нижнюю границу для шага h, в то время как для неявных методов (17) имеет место
верхняя граница, задаваемая условием сходимости итерационного процесса для решения нелинейных уравнений в неявных методах [6?:
h<
L maxi
1
P
j
|aij |
,
где L постоянная Липшица, которая задает неравенство
|f(t, x) ? f(t, y)| ? L|x ? y|
(35)
P
для любых x и y. Если положить maxi j |aij | = 1 (в действительности максимум
близок к единице), то получим следующее ограничение на шаг интегрирования: h < 1/L
или h? = n/L. Оценим постоянную Липшица L для исследуемой задачи. ассмотрим
отношение
|?f|2
n4 |?r|2 + |?v|2
=
,
|?x|2
|?r|2 + |?v|2
где ?x, ?r и ?v всевозможные разности векторов в соответствующих переменных.
Принимая |?r| = ? cos ? , |?v| = ? sin ? , где ? ? 0 и 0 ? ? ? ?/2, будем иметь
|?f|2
= (n4 ? 1) cos2 ? + 1.
|?x|2
Отсюда следует, что все значения отношения лежат между 1 и n4 . Следовательно, согласно (35) в качестве постоянной Липшица можно выбрать L = max(1, n2 ). Тогда
получим верхнюю границу шага
h? <
n
? 1.
max(1, n2 )
Таким образом, в лучшем случае, а именно при n = 1, когда верхняя граница максимальна, шаг интегрирования h? должен удовлетворять неравенствам
?(p + 1)! < hp+1
< 1.
?
Очевидно, условие
?(p + 1)! > 1
(36)
В. А. Авдюшев
40
означает, что порядок метода завышен и использование такого метода при вычислениях
в ариметике с точностью ? не логично в том смысле, что ту же точность результатов интегрирования можно получить с использованием методов более низких порядков.
Оптимальные порядки p неявных методов унгеКутты для различных ?, соответствующих одинарной, двойной, расширенной и четверной точности, представлены в табличном виде (следует иметь в виду, что эти порядки получены для задачи с n = 1, в ином
случае они могут быть меньше):
?
p
1.1 · 10?7
9
2.2 · 10?16
16
1.1 · 10?19
19
1.9 · 10?34
29
7. Фортран-код интегратора ауссаЭверхарта
Представленный выше ортран-код интегратора ауссаЭверхарта GAUSS_15 (до 15
порядка), основанного на гауссовых разбиениях Лобатто и адо, а также Лежандра,
доступен в интернете по адресу http://astro.tsu.ru/sh/Gauss_15.txt.
GAUSS_15(X,TS,TF,STEP,ERR,N,NOR,NI,NS,NBS,NF,FUN) заголовок процедуры интегрирования. Опишем входные и выходные параметры интегратора, помечая их соответственно индексами I и O . Итак, XIO интегрируемые переменные (до выполнения
процедуры начальные для значения TS, после конечные для TF); TSI и TFI начальное и конечное значения независимой переменной соответственно; STEPIO величина стартового шага интегрирования (если STEP = 0, то стартовый шаг выбирается
по ормуле (29)), после выполнения процедуры STEP величина предпоследнего шага;
ERRI значение ketol k для выбора переменного шага (27) (условие ERR 6= 0 задает режим переменного, ERR = 0 постоянного шага); NI число уравнений; NORI порядок
интегратора p; NII максимальное число итераций на шаге (при NI ? 0 итерационный
процесс выполняется до сходимости); NSO число шагов интегрирования за выполнение процедуры; NBSO число шагов, на которых итерационный процесс не сходится
(несходимость итерационного процесса регистрируется, когда число итераций достигает
100); NFO число обращений к процедуре правых частей; FUNI название процедуры
правых частей f .
8. Тестирование интегратора ауссаЭверхарта
в задаче двух тел
Интегратор тестировался на диеренциальных уравнениях задачи двух тел (30). Начальные условия задачи
r1 = 1 ? e,
r2 = 0,
v1 = 0,
v2 =
p
(1 + e)/(1 ? e)
соответствуют однопараметрическому семейству орбит с эксцентриситетом e в качестве
параметра и с единичной большой полуосью a. Интегрирование выполнялось на интервале 1000 оборотов для различных эксцентриситетов. При этом исследовались те или
иные характеристики интегратора.
Интегратор ауссаЭверхарта
8.1. Круговой случай
41
e=0
Как уже было замечено, шаг в интеграторе выбирается таким образом, чтобы сохранялась величина члена порядка k + 1. Если бы интегратор имел порядок k , то таким
способом можно было бы обеспечить сохранение локальной точности на всем интервале
интегрирования. Однако порядок интегратора ауссаЭверхарта существенно выше k .
На примере круговой задачи (e = 0) исследовалась зависимость реальной локальной
точности от задаваемой для выбора переменного шага.
В (16) величину k + 1 члена kek k можно оценить как
kek k ?
h
hk+1
kAk k ?
kx(k+1) k,
k+1
(k + 1)!
(37)
тогда как ошибка численного решения порядка p будет
kep k ?
hp+1
kx(p+1) k.
(p + 1)!
(38)
Поскольку решение x имеет вид (31), то
?
kx(i) k = a 1 + n2 ni .
(39)
Подставляя (39) в (37) и (38), получим
?
p?k
p+1
(a 1 + n2 )? k+1
kep k ?
((k + 1)!kek k) k+1 .
(p + 1)!
В частности, для разбиения ауссаадо (p = 2k + 1)
k+1
Y
1
i
kep k ? ?
kek k2 .
2
a 1 + n i=1 k + 1 + i
(40)
0.001
????????? ???????? |'x|
1E-005
1E-007
1E-009
1E-011
1E-013
1E-015
1E-017
1E-010
1E-008
1E-006
0.0001
0.01
ERR
ис. 1. Зависимость локальной точности от параметра ERR
1
В. А. Авдюшев
42
На рис. 1 (см. с. 41) представлены оценки реальной локальной ошибки
kep k = k?xk =
q
?r12 + ?r22 + ?v12 + ?v22
в зависимости от ERR (? kek k) для p = 11. Оценивание выполнялось в режиме
STEP = 0 на первом шаге, подбираемом в соответствии с параметром ERR. Видно, что
полученные результаты хорошо согласуются с (40). В частности, экспериментальная
оценка при ERR > 10?6 (ERR < 10?6 область доминирования вычислительных ошибок) подтверждает квадратичную зависимость kep k ? kek k2 , отчего главным образом
реальная точность численного решения существенно выше задаваемого значения параметра ERR.
8.2. Слабоэксцентричный случай
e = 0.1
Для оценки оптимизации составленного ортран-кода сравнили быстродействия (в затратах процессорного времени) интегратора GAUSS_15 и широко используемого в небесной механике интегратора RADAU_27. При этом в обоих случаях выполнялись все
операции, определяемые основными ормулами. В то же время расхождение двух численных решений, полученных интеграторами, оказалось меньше методической ошибки
на несколько порядков.
На рис. 2 показаны отношения временных затрат, которые потребовались интеграторам для получения численных решений 7, 11 и 15 порядков. Как видно, интегратор
GAUSS_15 работает быстрее, однако этот выигрыш с увеличением порядка уменьшается. Впрочем, следует заметить, что в задачах со сложной ункцией f оперативность
интеграторов должна быть одинаковой, поскольку большая часть времени будет затрачиваться на процедуру ее вычисления FUN и тогда быстродействие главным образом
определяется числом обращения к этой процедуре.
??????? ?? ???????? tRadau_27/tGauss_15
1.3
1.2
1.1
1
7
11
??????? ?????? p
15
ис. 2. Оптимальность интегратора GAUSS_15 относительно RADAU_27
Интегратор ауссаЭверхарта
43
p
Далее была оценена ошибка интегрирования |?r| = ?r12 + ?r22 в зависимости от
величины постоянного шага h для порядков 211. Интегрирование выполнялось с шестью итерациями (на шаге) на интервале 1000 оборотов.
В случае k = 13 представленные на рис. 3 характеристики показывают, что при достаточно больших h (глобальная) ошибка решения p-порядка, как и ожидалось, с уменьшением шага ведет себя как |?r| ? hp , что указывает на соответствие ошибки методу
порядка p. Для k = 4 и 5 ситуация несколько иная. Это объясняется тем, что шести
итераций недостаточно для сходимости итерационного процесса на шаге и получаемое
решение не соответствует порядку интегрирующих ормул. Случайное поведение характеристик ниже |?r| = 10?6 обусловлено влиянием ошибок округления.
Следует заметить, что несмотря на крутой скат характеристики нечетного порядка
2k + 1 в сравнении с характеристикой четного порядка 2k (для k = 13) определенная
точность для четного порядка достигается при большем шаге, т. е. быстрее. Это связано
с тем, что при симметричном разбиении ауссаЛобатто интегратор ауссаЭверхарта
обладает геометрическими свойствами и методическая ошибка вдоль независимой переменной t эволюционирует медленнее, чем при разбиении ауссаадо.
На рис. 4 показано поведение ошибки с изменением t для четных и нечетных порядков (k = 35). Приведенные результаты были получены в режиме NI ? 0, т. е. до
полной сходимости итерационного процесса. При этом на шаге требовалось от 11 до 16
итераций. Как следует из рисунка, при постоянном шаге (h = 2?/16) ошибка интегрирования для четных порядков ведет себя линейно, тогда как для нечетных порядков квадратически. Даже несмотря на то что интегратор порядка 2k+1 в начале интегрирования дает решение точнее интегратора порядка 2k , в конце интегрирования вследствие
разного роста ошибок первый уже уступает второму в точности почти на два порядка.
Очевидно, с увеличением интервала интегрирования это преимущество для каждого интегратора с симметричным разбиением будет возрастать. К сожалению, такие свойства
интеграторов четных порядков имеют место только при использовании постоянного
0.1
3
2
0.01
???????? |'r|
0.001
5
4
11
0.0001
1E-005
7
6
10
8
1E-006
1E-007
9
1E-008
1E-009
0.001
0.01
0.1
??? ?????????????? h
ис. 3. Зависимость точности от величины шага интегрирования
1
В. А. Авдюшев
44
0.1
0.01
0.001
???????? |'r|
0.0001
1E-005 6
1E-006
7
8
1E-007
1E-008
1E-009
??????? ??????
???????? (h const)
?????? (h const)
?????? (h zconst)
9
1E-010 10
1E-011
1E-012 11
1E-013
1
10
??????
100
1000
ис. 4. Поведение ошибки для решений четных и нечетных порядков
ис. 5. Характеристики точностьбыстродействие в зависимости от числа итераций на шаге
шага [5?. При автоматическом выборе шага ошибка так же, как и для интеграторов
нечетных порядков, ведет себя квадратичным образом.
Далее был исследован вопрос о соотношении точности и быстродействия при увеличении числа итераций на шаге. Путем вариации ERR оценивались точность и быстродействие интегратора 11 порядка с переменным шагом на интервале 1000 оборотов для
NI = 16. езультаты приведены на рис. 5. Видно, что уже на второй итерации можно
получить довольно хорошее решение, хотя соответствующая характеристика заметно
отклоняется от степенной зависимости. Поэтому для уверенного результата необходимо использовать, по крайней мере, три итерации на шаге. Следует также заметить, что
увеличение числа итераций не существенно понижает эективность интегрирования.
Интегратор ауссаЭверхарта
45
8.3. Сильноэксцентричный случай
e = 0.9
Очевидно, сильноэксцентричные орбиты необходимо интегрировать с переменным шагом. Поведение выбираемого шага на одном обороте для e = 0.9 представлено на рис. 6,
где также приведены ункции, пропорциональные |r|, |r|3/2 , |r|2 и |v|?1 . Как видно,
изменение шага очень близко к поведению |r|3/2 . Это означает, что шаг переменной величины по t соответствует почти постоянному шагу по так называемой эллиптической
аномалии.
1
??? ?????????????? h
~ |r |2
0.1
|r|3/2 ~
~ |r|
0.01
~ |v|1
0.001
0
40
80
120
????? ????
160
200
ис. 6. Выбор шага интегрирования в случае сильновытянутой орбиты
1
???????? |'r|
0.1
??????????
GAUSS_15
RADAU_27
0.01
0.001
0.0001
1E-005
1E-006
1000000
10000000
?????????????? NF
100000000
ис. 7. Характеристики точностьбыстродействие для интеграторов GAUSS_15 и RADAU_27
в экстраэксцентричном случае
В. А. Авдюшев
46
8.4. Экстраэксцентричный случай
e = 0.999
Наконец, был исследован алгоритм выбора шага при долгосрочном интегрировании
очень вытянутой орбиты с эксцентриситетом e = 0.999. Интегрирование выполнялось
на интервале 1000 оборотов двумя интеграторами: GAUSS_15 и RADAU_27. езультаты показали (рис. 7), что RADAU_27 заметно менее эективен, нежели GAUSS_15.
Так, при одинаковом быстродействии первый интегратор дает решение с точностью ниже на порядок и более. Это связано с тем, что в RADAU_27 заложен неверный алгоритм
выбора шага для интегрирования систем с уравнениями первого порядка, а именно, алгоритм не точно соответствует ормуле (27): вместо степени 1/(k + 1) в RADAU_27
используется 1/(k + 2). Очевидно, эта ошибка легко устраняется.
Таким образом, в работе представлена новая версия интегратора ауссаЭверхарта
GAUSS_15. На примере плоской задачи двух тел при различных начальных условиях
экспериментально показано, что предложенная программная реализация интегратора
значительно эективнее ранее используемой, а также предоставляет больше возможностей для высокоточного и оперативного численного интегрирования.
Список литературы
[1?
[2?
Everhart E.
P. 389.
A new method for integrating orbits // Bull. Amer. Astronom. So. 1973. Vol. 5.
Everhart E. Impliit single sequene methods for integrating orbits // Celest. Meh. 1974.
Vol. 10. P. 3555.
[3?
Everhart E.
[4?
Buther J.C.
[5?
[6?
An eient integrator that uses GaussRadau spaings // Dynamis of Comets:
Their Origin and Evolution. Pro. of IAU Colloq. 83. Italy, 1984 / Eds. A. Carusi and
G.B. Valsehi. Dordreht: Reidel, Astrophys. and Spae Siene Library Rame, 1985. Vol. 115.
P. 185202.
Impliit RungeKutta proesses // Math. Comput. 1964. Vol. 18. P. 5064.
Hairer E., Lubih C., Wanner G. Geometri Numerial Integration: Struture-Preserving
Algorithms for Ordinary Dierential Equations / Springer Series in Comput. Math. Springer,
2002. 536 p.
Solving Ordinary Dierential Equations. Nonsti
Problems / Springer Series in Comput. Math. Springer, 1993. 544 p.
Hairer E., Norsett S.P., Wanner G.
Поступила в редакцию 27 августа 2009 г.,
с доработки 26 сентября 2009 г.
Документ
Категория
Без категории
Просмотров
25
Размер файла
496 Кб
Теги
интеграторах, эверхарта, гаусса
1/--страниц
Пожаловаться на содержимое документа