close

Вход

Забыли?

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

?

Алгоритм дифференцирования основанный на методе дополнительных переменных.

код для вставкиСкачать
УДК 517.9:519.6
Вестник СПбГУ. Сер. 10, 2013, вып. 2
К. М. Брэгман
АЛГОРИТМ ДИФФЕРЕНЦИРОВАНИЯ, ОСНОВАННЫЙ
НА МЕТОДЕ ДОПОЛНИТЕЛЬНЫХ ПЕРЕМЕННЫХ
Введение. Средства символьного и автоматического дифференцирования функций
многих переменных [1, 2] настолько важны и востребованы, что входят во все основные
математические пакеты компьютерной алгебры (см., например, [3]). В настоящей статье
предлагается новый алгоритм нахождения производных системы функций многих переменных весьма общего вида. Речь идет о функциях, которые можно получить из своих
аргументов при помощи конечного числа операций +, −, ×, / и суперпозиций функций,
удовлетворяющих полным полиномиальным системам уравнений в частных производных (и, в частности, системам обыкновенных дифференциальных уравнений – ОДУ).
Класс таких функций включает в себя элементарные функции, специальные функции
математической физики (за исключением теоретико-числовых функций) и многие другие. Алгоритм основан на методе дополнительных переменных (МДП) и использует
библиотеку, содержащую набор таких функций (а точнее, полиномиальных дифференциальных уравнений для этих функций). Как и в статьях [4, 5], предполагается, что
исходная система функций и функции библиотеки могут зависеть также и от конечного числа параметров. Последовательное введение новых (дополнительных) переменных сводит все функции исходной системы к полиномиальным функциям по исходным
и новым аргументам, причем новые аргументы удовлетворяют полной полиномиальной системе уравнений в частных производных (по исходным аргументам). Получение
производных исходной системы сводится после этого к использованию «цепного правила» дифференцирования сложных функций и, в итоге, к перемножению и сложению
полиномов, а результаты выражаются полиномами по исходным и новым переменным.
Для удобства читателя в начале статьи (в п. 1) сразу приводится нетривиальный пример, поясняющий предлагаемый алгоритм и его отличие от известных алгоритмов символьного и автоматического дифференцирования: рассматриваются две функции трех
переменных и пяти параметров и показывается, как после введения нескольких дополнительных переменных вычисление частных производных любого порядка для этих
функций сводится к сложению и умножению полиномов. В п. 2 приводятся необходимые сведения о классах функций, дифференциальных уравнениях и библиотеках
функций, в п. 3 описывается сам алгоритм, а в п. 4 – пример его применения.
Отметим наконец, что для нахождения производных функций многих переменных
мы применяем те же средства, что и в статье [5], поэтому здесь, насколько возможно, используем аналогичные обозначения, а также аналогичные последовательность
и способ изложения. Это позволило предоставить пользователю возможности преобразования дифференциальных уравнений к полиномиальной форме и дифференцирования функций в рамках схожих алгоритмов метода дополнительных переменных. Естественно, что читателю рекомендуется сначала ознакомиться со статьей [5] (хотя это
и не обязательно), а затем уже с настоящей работой, содержащей алгоритм еще одного
применения МДП.
Брэгман Константин Михайлович – аспирант, Санкт-Петербургский государственный университет; e-mail: kostjab@yandex.ru.
c К. М. Брэгман, 2013
14
1. Поясняющий пример. Рассмотрим систему из двух функций трех переменных
x cos(ωx2 )
y0,1 = x1 3
(ax1 + bx2 )7 x3 sin(ωx2 ),
y0,2 = ln12 x1 (c3 x1 x2 x3 cos3 (ωx2 ) + coth(cd d) sin(ωx2 ))3 ,
где x1 , x2 , x3 – переменные, а a, b, c, d, ω – параметры. Вводя еще переменные
ϕ1 = sin(ωx2 ),
ϕ2 = cos(ωx2 ),
ϕ3 = ln x1 ,
ϕ4 = x−1
1 ,
x cos(ωx2 )
ϕ5 = x1 3
,
вместо y0,1 , y0,2 получаем полиномы (по переменным x = (x1 , x2 , x3 ), ϕ = (ϕ1 , . . . , ϕ5 ))
y1 = (ax1 + bx2 )7 x3 ϕ1 ϕ5 ,
3
3
d
3
y2 = ϕ12
3 (c x1 x2 x3 ϕ2 + coth(c d)ϕ1 ) ,
причем ϕ является решением полной системы:
∂ϕ1 /∂x1 = 0,
∂ϕ1 /∂x2 = ωϕ2 ,
∂ϕ2 /∂x2 = −ωϕ1 ,
∂ϕ3 /∂x3 = 0,
∂ϕ2 /∂x3 = 0,
∂ϕ4 /∂x1 = −ϕ24 ,
∂ϕ5 /∂x1 = x3 ϕ2 ϕ4 ϕ5 ,
∂ϕ1 /∂x3 = 0,
∂ϕ2 /∂x1 = 0,
∂ϕ3 /∂x1 = ϕ4 ,
∂ϕ4 /∂x2 = 0,
∂ϕ5 /∂x2 = −ωx3 ϕ1 ϕ3 ϕ5 ,
∂ϕ3 /∂x2 = 0,
∂ϕ4 /∂x3 = 0,
∂ϕ5 /∂x3 = ϕ2 ϕ3 ϕ5 .
Если использовать обозначения
Di y0,r = ∂ i1 +i2 +i3 y0,r /∂xi11 ∂xi22 ∂xi33 ,
e1 = (1, 0, 0),
i = (i1 , i2 , i3 ),
e2 = (0, 1, 0),
ik ∈ [0 : +∞),
e3 = (0, 0, 1),
то для вычисления искомых производных можно использовать рекуррентные формулы
(«цепное правило» дифференцирования сложных функций)
Dek y0,r = ∂yr /∂xk +
5
∂yr /∂ϕm · ∂ϕm /∂xk ,
m=1
D
i+ek
i
y0,r = ∂D y0,r /∂xk +
5
∂Di y0,r /∂ϕm · ∂ϕm /∂xk ,
r = 1, 2;
k = 1, 2, 3,
m=1
причем каждая из получаемых величин Dek y0,r , . . . , Di+ek y0,r будет полиномом
по x1 , x2 , x3 , ϕ1 , . . . , ϕ5 , коэффициенты которого зависят от α = (a, b, c, d, ω).
Для вычисления производной ∂ i1 +i2 +i3 y0,r /∂xi11 ∂xi22 ∂xi33 для любых значений x, α
можно сначала определить ϕ, а затем значение полинома Di y0,r , причем его можно
рассчитать по явной формуле (если она получена заранее по приведенным выше рекуррентным соотношениям) или непосредственно по рекуррентным соотношениям.
Таким образом, можно сказать, что предлагаемый алгоритм использует, как и общепринятые методы символьного и автоматического дифференцирования, «цепное правило», но благодаря технике МДП (величины ∂yr /∂xk , ∂yr /∂ϕm , ∂ϕm /∂xk в этом случае
оказываются полиномами по исходным и дополнительным переменным) алгоритм применения этого правила сводится сначала к вычислению дополнительных переменных
и затем к перемножению и сложению полиномов.
15
Заметим, что алгоритм не отвечает на два ключевых вопроса: как автоматизировать
выбор дополнительных переменных и как найти для них полную систему уравнений
в частных производных с полиномиальными правыми частями. Ответы на эти вопросы
будут даны в п. 3.
2. Дифференциальные уравнения, классы функций, библиотеки. Здесь
приводятся необходимые сведения из статей [4, 5]. В п. 2.1 вводятся в рассмотрение
полные полиномиальные системы дифференциальных уравнений и классы функций,
которые таким системам удовлетворяют; в п. 2.2 – классы функций многих переменных, на которые ориентирован предлагаемый в настоящей статье алгоритм (см. п. 3).
В п. 2.3 обсуждается основной используемый в этом алгоритме инструмент – библиотека.
2.1. Дифференциальные уравнения, классы Σσ . Положим
x = (x1 , . . . , xm ) ∈ Rm ,
t = (t1 , . . . , ts ) ∈ Rs ,
α = (α1 , . . . , αω ) ∈ Rω
и рассмотрим систему уравнений
∂xi /∂tj = fij (x, α),
i ∈ [1 : m],
j ∈ [1 : s],
где все функции fij – полиномы по x1 , . . . , xm с коэффициентами, зависящими от α.
Такую систему называют полиномиальной. В случае s = 1 это система ОДУ.
Говорят, что скалярная функция ϕ аргумента x = (x1 , . . . , xσ ) удовлетворяет полиномиальной системе, если она является одной из компонент вектор-функции (того
же аргумента x = (x1 , . . . , xσ )) – решения некоторой полиномиальной системы. Класс
скалярных функций аргумента x = (x1 , . . . , xσ ), удовлетворяющих полиномиальной системе, обозначим Σσ . Ясно, что Σ1 ⊂ Σ2 ⊂ . . . . Все элементарные функции и большое
число специальных функций принадлежат Σ1 (и тем более Σσ при σ > 1).
σ
σ
2.2. Классы Fm
. Классом Fm
, σ, m ∈ [1 : +∞), называют множество скалярных
функций аргумента x = (x1 , . . . , xm ), которые можно получить из x1 , . . . , xm при помощи конечного числа операций +, −, ×, / и конечного числа функций ϕ1 , . . . , ϕl ∈ Σσ
и их суперпозиций.
2.3. Расширения и библиотеки. Если ϕ ∈ Σσ , то существует вектор-функция
(ϕ1 , . . . , ϕn ) – решение полиномиальной системы, где ϕ1 = ϕ. Расширением ϕ называют
множество {ϕ1 , . . . , ϕn }. Объединение расширений нескольких функций называют библиотекой. Фактически, библиотека содержит набор полиномиальных систем, а не задач Коши. Всякое подмножество библиотеки называют подбиблиотекой, если оно само
является библиотекой. Библиотеку называют автономной или неавтономной соответственно тому, какими уравнениями она представлена – автономными или неавтономными. Естественно пользоваться как автономными, так и неавтономными библиотеками,
хотя можно было бы ограничиться автономными библиотеками, так как неавтономную
библиотеку всегда можно свести к автономной.
Разделом называют подбиблиотеку, которая не пересекается ни с одной другой
подбиблиотекой, не являющейся ее частью. Объединение разделов является разделом. Раздел называют простым, если он не содержит других разделов. Библиотека
является объединением своих разделов, она делится на разделы и сама может быть
своим единственным разделом. Для разбиения библиотеки на разделы и их нумерации можно предложить различные алгоритмы, основанные на очевидном утверждении: если в списке расширения функции есть функция из некоторого раздела, то и сама
функция принадлежит этому разделу. При пополнении библиотеки новыми функциями
16
имевшуюся нумерацию разделов можно сохранить, а новым разделам (если появятся)
присвоить новые номера (старые разделы могут пополниться новыми функциями).
Для пользователя библиотека – пополняемая таблица, состоящая из потенциально
неограниченного количества строк и одиннадцати столбцов, полное описание которых
дано в п. 2.1 статьи [5].
3. Алгоритм дифференцирования. Рассмотрим систему функций
yr = fr (x),
r ∈ [1 : N ] (или y = f (x)),
(1)
при x = (x1 , . . . , xm ) ∈ Rm , y = (y1 , . . . , yN ), f = (f1 , . . . , fN ) ∈ RN и при условии, что
σ
все fr (x) принадлежат классу Fm
и записаны в терминах функций имеющейся в нашем
распоряжении библиотеки. Кроме того, будем считать, что fr (x) могут зависеть от параметров α = (α1 , . . . , αω ) ∈ Rω , хотя мы и не указываем на это явно в обозначениях.
Мы собираемся здесь вывести формулы для производных ∂ i1 +···+im yr /∂xi11 . . . ∂ximm .
В п. 3.1 покажем, как, вводя дополнительные переменные при помощи библиотеки,
преобразовать систему к полиномиальной форме, и выведем рекуррентные формулы
для производных введенных переменных по x1 , . . . , xm . В п. 3.2 получим формулы для
искомых производных, а в п. 3.3 опишем схему алгоритма.
3.1. Элементарное преобразование системы. Преобразование исходной системы к полиномиальной состоит из конечного числа K последовательных шагов, причем на каждом из них проводится однотипное «элементарное» преобразование, которое мы далее рассмотрим. Ради симметрии обозначений, для величин yr , y, fr , f, x, m
в (1) будем использовать также и обозначения y0,r , y0 , f0,r , f0 , x0 , m(0) соответственно.
На каждом шаге вводятся дополнительные переменные, причем в новых переменных
система (1) представляется в виде (k – номер шага)
yk,r = fk,r (xk ),
r ∈ [1 : N ] или y = fk (xk ),
k ∈ [1 : K],
(2)
xk = (xk−1 ; X k ) = (x1 , . . . , xm(k−1) ; xm(k−1)+1 , . . . , xm(k−1)+n(k) ) = (x1 , . . . , xm(k) ),
m(k) = m(k − 1) + n(k),
0
x = x = (x1 , . . . , xm(0) ),
m(0) = m,
k ∈ [1 : K],
y0,r = yr ,
f0,r = fr ,
а натуральные числа n(k) и функции f1,r , f2,r , . . . определяются на каждом очередном
шаге. При k = 0 системы (1) и (2) совпадают. Первый шаг состоит в переходе от системы (1) к системе (2) при k = 1, а на последнем шаге, при k = K, система (2) становится
полиномиальной. Перейдем к описанию k-го шага (k ∈ [1 : K]).
а). Ищем в выражениях fk−1,r (xk−1 ), r ∈ [1 : N ] функцию вида ψk (P k (xk−1 ))ρ(k) ,
где P k (xk−1 ) = (Pk,1 (xk−1 ), . . . , Pk,ϑ(k) (xk−1 )), Pk,ν – алгебраические полиномы, ψk =
ψk (p1 , . . . , pϑ(k) ) – функция из библиотеки, а ρ(k), ϑ(k) ∈ [1 : +∞).
б). Пусть ψk = ϕ1 , а ϕ1 , . . . , ϕn(k) – все функции расширения функции ψ k и пусть
∂ϕl /∂pν = Qνk,l (pk ; ϕk ),
l ∈ [1 : n(k)],
ν ∈ [1 : ϑ(k)],
(3)
где pk = (p1 , . . . , pϑ(k) ), ϕk = (ϕ1 , . . . , ϕn(k) ), а Qνk,l (pk ; ϕk ) – полиномы по всем n(k) +
ϑ(k) своим аргументам. Вводим n(k) новых дополнительных переменных
xm(k−1)+1 = ϕ1 (P k (xk−1 )), . . . , xm(k) = ϕn(k) (P k (xk−1 )).
(4)
17
в). Заменяя в fk−1, r (xk−1 ), r ∈ [1 : N ], все функции ϕ1 (P k (xk−1 )), . . . , ϕn(k) (P k (xk−1 ))
на новые переменные xm(k−1)+1 , . . . , xm(k) соответственно, получаем
r ∈ [1 : N ].
yk,r = fk,r (xk ),
(5)
г). При l ∈ [1 : n(k)], s ∈ [1 : m], выписываем производные (см. (3), (4))
ϑ(k)
d
d
xm(k−1)+l (xk−1 ) =
Qνk, l (P k (xk−1 ); X k ) ·
Pk,ν (xk−1 ),
dxs
dx
s
ν=1
(6)
m(k−1)
∂
d
dxη (xk−1 )
k−1
Pk,ν (x
)=
Pk,ν (xk−1 ) ·
,
dxs
∂xη
dxs
η=1
(7)
причем d/dxs – производная сложной функции по аргументу xs , s ∈ [1 : m] при постоянных значениях аргументов xj , j ∈ [1 : m], j = s, а xη (xk−1 ) = xη при η ∈ [1 : m].
Для того чтобы воспользоваться формулами (6), (7) для последовательного (при
k = 1, . . . , K) нахождения производных дополнительных переменных по аргументам x1 , . . . , xm , необходимо на каждом шаге иметь в своем распоряжении величины
ϑ(k),
n(k), Qνk,l (p1 , . . . , pϑ(k) ; ϕ1 , . . . , ϕn(k) ),
Pk,ν (xk−1 ),
l ∈ [1 : n(k)], ν ∈ [1 : ϑ(k)],
(8)
которые уже известны (см. a), б)).
3.2. Производные исходной системы. Если использовать обозначения
Di yr = ∂ i1 +···+im fr (x)/∂xi11 . . . ∂ximm ,
i = (i1 , . . . , im ),
ij ∈ [0 : +∞),
e1 = (1, 0, . . . , 0), . . . , em = (0, . . . , 0, 1),
то, применяя «цепное правило» дифференцирования, получаем
Dej yr =
Di+ej yr =
∂
fK,r (xK ) +
∂xj
∂Di yr
+
∂xj
m(K)
μ=m(0)+1
m(K)
μ=m(0)+1
∂Di yr dxμ
·
,
∂xμ
dxj
∂fK,r (xK ) dxμ
·
,
∂xμ
dxj
r ∈ [1 : N ],
j ∈ [1 : m],
причем производные dxμ /dxj вычисляются по формулам (6), (7) и dxμ /dxj = δμ,j при
μ, j ∈ [1 : m].
Найдем более удобный вариант этих формул, для чего полиномы yr в (5) при k = K
запишем в виде
u
yr =
ar,η xλ(η) , x = (x1 , ..., xm (K) ),
(9)
η=0
λ(η) = (λη,1 , . . . , λη,m (K) ), r ∈ [1 : N ],
где λ(0) = (0, . . . , 0), xλ(1) = x1 , . . . , xλ(m(K)) = xm(K) , а xλ(m(K)+1) , . . . , xλ(u) – все
различные нелинейные мономы, из которых составлены полиномы fK,r , r ∈ [1 : N ].
18
Полагая D(0,...,0) yr = yr , обобщим представление (9) на Di yr :
i
D yr =
u(i)
ar,η (i)xλ(i,η) ,
x = (x1 , ..., xm(K) ),
η=0
λ(i, η) = (λη,1 (i), . . . , λη,m (K) (i)),
r ∈ [1 : N ].
Обозначения здесь аналогичны и очевидны, причем эта формула включает в себя (9)
как частный случай при i = (0, ..., 0). Дифференцируя ее, имеем следующие рекуррентные соотношения для нахождения искомых производных:
⎞
⎛
u(i)
m(K)
dx
μ
⎠.
Di+ej yr =
ar,η (i) ⎝
λη,μ (i)xλ(i, η)−eμ ·
(10)
dxj
η=0
μ=1
Вычисления по формулам (10) могут отличаться тем, в какой последовательности находятся производные (полиномы по исходным и введенным аргументам). Производные
dxμ /dxj находятся заранее по рекуррентным формулам (6), (7).
3.3. Схема алгоритма.
Шаг 1. Заносим систему (1) в таблицу (например, в пакете Wolfram Mathematica).
Введенная система может содержать параметры. Если для обозначения каких-то
из этих параметров использованы символы вида αγ , где γ – натуральное число, то полагаем, что Ω равно максимальному из таких γ, а если таких параметров нет, то полагаем Ω = 0. Параметры, от которых зависят функции подбиблиотеки, переобозначим
символами αΩ+ 1 , αΩ+ 2 , . . . и запомним таблицу соответствий новых и старых обозначений. Образуем подбиблиотеку, состоящую из объединения расширений всех функций,
участвующих в написании исходной системы.
Шаг 2. Преобразуем систему шаг за шагом по алгоритму п. 3.1, получая на каждом
из них систему вида (5), пока она не окажется полиномиальной.
Шаг 3. Вычисляем производные по формулам (6), (7).
Заметим, что в программе можно эти производные вычислять и последовательно
в рамках шага 2.
Шаг 4. Расcчитываем производные исходной системы по алгоритму п. 3.2.
Приведенная на рисунке блок-схема отражает описанный выше алгоритм получения
производных функций системы (1).
4. Пример. Чтобы ограничиться той же библиотекой, что и в статье [5], рассмотрим
систему из шести функций, которые немногим отличаются от правых частей полной
системы уравнений, использованной там в качестве примера:
y1 = f1 (x1 , x2 , x3 ) = x33 (sin cos(a ln2 x2 + bx3 ) + b ln4 x2 ) +
+ sin5 (a ln2 x2 + bx3 ) + Dv(g 2 , 1, −1; x1 ),
y2 = f2 (x1 , x2 , x3 ) = (sin cos(a ln2 x2 + bx3 ) + b ln4 x2 ) +
+ sin5 (a ln2 x2 + bx3 ) + Dv(g 2 , 1, −1; x1 ),
y3 = f3 (x1 , x2 , x3 ) = x22 cos sin(a ln2 x2 + bx3 ) +
+ cos4 (a ln2 x2 + bx3 ) + EK 2 (sin sin(a ln2 x2 + bx3 ), x1 ),
y4 = f4 (x1 , x2 , x3 ) = x22 cos sin(a ln2 x2 + bx3 ) +
+ cos4 (a ln2 x2 + bx3 ) + EK(sin sin(a ln2 x2 + bx3 ), x1 ),
(11)
19
Сообщение
об ошибке
Вход
и заведение системы
6
Нет
?
Заведение
и проверка
данных
Возможно ли построение подбиблиотеки
Да
?
Построение подбиблиотеки для исходной системы.
Формирование разделов и таблицы параметров
Подготовка
?
Является ли полиномиальной полученная
система?
6
Нет -
Найти в системе
функцию, аргументы
которой – полиномы
?
Замена в системе функций
расширения новыми переменными
Преобразо?
Да
вания
Нахождение производных для дополнительных
переменных
?
Нахождение искомых производных исходной системы
?
Вывод исходной системы функций, таблицы параметров,
дополнительных переменных и их производных,
системы полиномов и их производных
Блок-схема алгоритма дифференцирования
20
Вывод
данных
y5 = f5 (x1 , x2 , x3 ) = x21 (ch2 (a ln2 x2 + bx3 ) +
+ sin(a ln2 x2 + bx3 )) + sh5 (a ln2 x2 + bx3 ),
y6 = f6 (x1 , x2 , x3 ) = x1 (ch2 (a ln2 x2 + bx3 ) +
+ sin(a ln2 x2 + bx3 )) + sh5 (a ln2 x2 + bx3 ),
где α = (α1 , α2 , α3 ) = (a, b, g) – вещественные постоянные (параметры).
Упомянутая выше библиотека является объединением шести подбиблиотек (Пб), которые мы сейчас рассмотрим вместе с соответствующими дифференциальными уравнениями. Функции обозначаются ϕ, ϕ1 , . . . , а их аргументы – p, p1 , . . . [5].
Пб1. Состоит из одной функции одного аргумента – ϕ = inv(p) = p−1 , которая
удовлетворяет уравнению dϕ/dp = −ϕ2 , а расширение этой функции состоит из нее
одной.
Пб2. Состоит из двух функций одного аргумента – ϕ1 = ln p, ϕ2 = inv(p), которые
удовлетворяют системе dϕ1 /dp = ϕ2 , dϕ2 /dp = −ϕ22 , причем расширение функции ϕ1
состоит из функций ϕ1 , ϕ2 , а расширение функции ϕ2 – из нее одной (т. е. совпадает
с Пб1).
Пб3. Состоит из двух функций одного аргумента – ϕ1 = sin p, ϕ2 = cos p, которые
удовлетворяют системе dϕ1 /dp = ϕ2 , dϕ2 /dp = −ϕ1 , а расширение каждой из функций
ϕ1 , ϕ2 состоит из функций ϕ1 , ϕ2 .
Пб4. Состоит из двух функций одного аргумента – ϕ1 = sh(p), ϕ2 = ch(p), которые
удовлетворяют системе dϕ1 /dp = ϕ2 , dϕ2 /dp = ϕ1 , а расширение каждой из функций
ϕ1 , ϕ2 состоит из функций ϕ1 , ϕ2 .
Пб5. Состоит из четырех функций двух аргументов – ϕ1 = EK(p1 , p2 ) = E(e, M ),
где EK = E – эксцентрическая аномалия, p1 = e – эксцентриситет, p2 = M – средняя
аномалия (см. [6]),
ϕ2 = EKs(p1 , p2 ) = sin ϕ1 ,
ϕ3 = EKc(p1 , p2 ) = cos ϕ1 ,
ϕ4 = EKi(p1 , p2 ) = (1 − p1 cos ϕ1 )−1 ,
которые удовлетворяют системе
∂ϕ1 /∂p1 = ϕ2 ϕ4 ,
∂ϕ1 /∂p2 = ϕ4 ,
∂ϕ2 /∂p1 = ϕ2 ϕ3 ϕ4 ,
∂ϕ3 /∂p1 =
−ϕ22 ϕ4 ,
∂ϕ2 /∂p2 = ϕ3 ϕ4 ,
∂ϕ3 /∂p2 = −ϕ2 ϕ4 ,
∂ϕ4 /∂p1 = ϕ3 ϕ24 − p1 ϕ22 ϕ34 ,
∂ϕ4 /∂p2 = −p1 ϕ34 ϕ2 ,
причем расширение функции ϕ1 состоит из ϕ1 , ϕ2 , ϕ3 , ϕ4 , а расширение каждой
из функций ϕ2 , ϕ3 , ϕ4 – из всех этих трех функций.
Пб6. Состоит из двух функций ϕ1 , ϕ2 аргумента p, причем первая из них – функция
Вебера Dv, а вторая – еe производная, которая обозначена DV . Как известно [7], они
удовлетворяют системе уравнений
dϕ1 /dp = ϕ2 ,
dϕ2 /dp = −(ap2 + bp + c)ϕ1 ,
где a, b, c – параметры. Расширение функций ϕ1 , ϕ2 состоит из функций ϕ1 , ϕ2 .
Переходим к реализации алгоритма, описанного в п. 3.
21
Шаг 1. Пользователь заводит систему (11) (например, в пакете «Mathematica»).
Шаг 2. Преобразуем систему шаг за шагом по алгоритму п. 3.1 (шаг 2.1, шаг 2.2,. . . ),
получая на каждом из них систему вида (5), пока она не окажется полиномиальной.
Всего шагов оказалось семь (т. е. K = 7). Ради экономии места, выпишем сразу все
дополнительные переменные, расположив их по строкам так, что первая из переменных
в k-й строке соответствует функции ψk = ψk (p1 , . . . , pϑ(k) ), найденной на k-м шаге, а все
переменные в этой строке соответствуют функциям ее расширения ψk = ϕ1 , . . . , ϕn(k)
(см. п. 3.1 а), б)):
x4 = ln x2 , x5 = x−1
2 ,
x14 = EK(x11 , x1 ),
x6 = sin(ax24 + bx3 ),
x7 = cos(ax24 + bx3 ),
x8 = sh(ax24 + bx3 ),
x9 = ch(ax24 + bx3 ),
x10 = cos x6 ,
x11 = sin x6 ,
x12 = sin x7 ,
x13 = cos x7 ,
x15 = EKs(x11 , x1 ),
2
x18 = Dv(g , 1, −1; x1 ),
x16 = EKc(x11 , x1 ),
(12)
x17 = EKi(x11 , x1 ),
2
x19 = DV (g , 1, −1; x1 ).
После седьмого шага исходная система (11) приобретает полиномиальную форму:
y7,1 = f7,1 (x1 , . . . , x19 ) = x33 (x12 + bx44 ) + x56 + x18 ,
y7,2 = f7,2 (x1 , . . . , x19 ) = (x12 + bx44 ) + x56 + x18 ,
y7,3 = f7,3 (x1 , . . . , x19 ) = x22 x10 + x47 + x214 ,
y7,4 = f7,4 (x1 , . . . , x19 ) = x22 x10 + x47 + x14 ,
y7,5 = f7,5 (x1 , . . . , x19 ) = x21 (x29 + x6 ) + x58 ,
y7,6 = f7,6 (x1 , . . . , x19 ) = x1 (x29 + x6 ) + x58 .
Шаги 3, 4. Для того чтобы воспользоваться формулами (6), (7) для последовательного (при k = 1, . . . , 7) нахождения производных x4 , . . . , x19 по x1 , x2 , x3 , необходимо
на каждом шаге иметь в своем распоряжении величины (8). Пользуясь библиотекой,
при пошаговом введении переменных (а они у нас так и записаны в (12)), их можно
получить также по шагам. Последовательность действий опишем на примере шестого
шага, на котором найдена функция ψ6 (P6,1 (x5 ), P6,2 (x5 ))2 = (EK(x11 , x1 ))2 .
Как видим, P6,1 (x5 ) = x11 , P6,2 (x5 ) = x1 , а ψ6 (p1 , p2 ) = EK(p1 , p2 ) – функция из подбиблиотеки Пб5. Так как расширение ψ6 (p1 , p2 ) = ϕ1 (p1 , p2 ), ϕ2 (p1 , p2 ) =
sin ϕ1 , ϕ3 (p1 , p2 ) = cos ϕ1 , ϕ4 (p1 , p2 ) = (1−p1 cos ϕ1 )−1 функции ψ6 (p1 , p2 ) состоит из четырех функций двух аргументов, то n(6) = 4, ϑ(6) = 2, а так как функции ϕ1 , . . . , ϕ4
удовлетворяют системе (Qi6,j = Qi6,j (p1 , p2 ; ϕ1 , . . . , ϕ4 ))
∂ϕ1
= Q16,1 = ϕ2 ϕ4 ,
∂p1
∂ϕ1
∂ϕ2
∂ϕ2
= Q26,1 = ϕ4 ,
= Q16,2 = ϕ2 ϕ3 ϕ4 ,
= Q26,2 = ϕ3 ϕ4 ,
∂p2
∂p1
∂p2
∂ϕ3
∂ϕ3
= Q16,3 = −ϕ22 ϕ4 ,
= Q26,3 = −ϕ2 ϕ4 ,
∂p1
∂p2
∂ϕ4
∂ϕ4
= Q16,4 = ϕ3 ϕ24 − p1 ϕ22 ϕ34 ,
= Q26,4 = −p1 ϕ34 ϕ2 ,
∂p1
∂p2
то и полиномы Qνk,l известны. Величины (8) для всех шагов получаются аналогично,
они представлены в табл. 1. В табл. 2 приведены используемые в формуле (7) частные
производные ∂Pk,ν (xk−1 )/∂xη , и, наконец, в табл. 3 даны рассчитанные по формулам
(6), (7) выражения (полиномы) для dxμ /dxj , μ = 1, . . . , 19; j = 1, 2, 3.
22
Теперь по формулам (10) можно получить выражения для производных функций
системы (11) по x1 , x2 , x3 . Пусть, например, требуется многократно вычислять производные ∂y6 /∂xj , j = 1, 2, 3. Выводя формулы для этих величин, имеем
∂y6 /∂x1 = (x29 + x6 ) + x1 (2x9 · dx9 /dx1 + dx6 /dx1 ) + 5x48 · dx8 /dx1 = x29 + x6 ,
∂y6 /∂x2 = dx1 /dx2 · (x29 + x6 ) + x1 (2x9 · dx9 /dx2 + dx6 /dx2 ) + 5x48 · dx8 /dx2 =
= 4ax1 x4 x5 x8 x9 + 2ax1 x4 x5 x7 + 10ax4 x5 x38 x9 = 2ax4 x5 (2x1 x8 x9 + x1 x7 + 5x38 x9 ),
∂y6 /∂x3 = dx1 /dx3 · (x29 + x6 ) + x1 (2x9 · dx9 /dx3 + dx6 /dx3 ) + 5x48 · dx8 /dx3 =
= x1 (2bx8 x9 + bx7 ) + 5bx48 x9 .
Таблица 1. Величины ϑ(k), n(k), Qνk,l (p1 , . . . , pϑ(k) ; ϕ1 , . . . , ϕn(k) ), Pk,ν (xk−1 )
k
ϑ (k)
n(k)
P k,ν
Qνk,l
1
1
2
P1,1 = x2
Q11,1 = ϕ2 ,
Q11,2 = −ϕ22
2
1
2
P2,1 = ax24 + bx3
Q12,1 = ϕ2 ,
Q12,2 = −ϕ1
3
1
2
P3,1 = ax24 + bx3
Q13,1 = ϕ2 ,
4
1
2
P4,1 = x6
Q14,1 = −ϕ2 ,
5
1
2
P5,1 = x7
Q15,1 = ϕ2 ,
Q13,2 = ϕ1
Q14,2 = ϕ1
Q15,2 = −ϕ1
6
2
4
P6,1 = x11 ,
P6,2 = x1
Q16,1 = ϕ2 ϕ4 ,
Q26,1 = ϕ4 ,
1
Q 6,2 = ϕ2 ϕ3 ϕ4 ,
Q26,2 = ϕ3 ϕ4 ,
Q16,3 = −ϕ22 ϕ4 ,
Q26,3 = −ϕ2 ϕ4 ,
Q16,4 = ϕ3 ϕ24 − p1 ϕ22 ϕ34 ,
Q26,4 = −p1 ϕ34 ϕ2
7
1
2
P7,1 = x1
Q17,1 = ϕ2 ,
Q17,2 = (1 − g 2 p2 − p)ϕ1
Таблица 2. Ненулевые величины ∂Pk,ν (xk−1 )/∂xη
k
1
2
3
4
5
6
7
∂ Pk,ν (xk−1 )/∂xη
∂P1,1 /∂x2 = 1
∂P2,1 /∂x3 = b, ∂P2,1 /∂x4 = 2ax4
∂P3,1 /∂x3 = b, ∂P3,1 /∂x4 = 2ax4
∂P4,1 /∂x6 = 1
∂P5,1 /∂x7 = 1
∂P6,1 /∂x11 = 1, ∂P6,2 /∂x1 = 1
∂P7,1 /∂x1 = 1
23
Таблица 3. Производные dxμ /dxj , μ = 1, . . . , 19; j = 1, 2, 3
μ
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
1
1
0
0
0
0
0
0
0
0
0
0
0
0
x17
x16 x17
−x15 x17
17
−x11 x15 x317
18
x19
19
(1 − g 2 x21 −
− x1 )x18
j
2
0
1
0
x5
−x25
2ax4 x5 x7
−2ax4 x5 x6
2ax4 x5 x9
2ax4 x5 x8
−2ax4 x5 x7 x11
2ax4 x5 x7 x10
−2ax4 x5 x6 x13
2ax4 x5 x6 x12
2ax4 x5 x7 x10 x17 (2x15 + 1)
2ax4 x5 x7 x10 x16 x17 (2x15 + 1)
−2ax4 x5 x7 x10 x15 x17 (2x15 + 1)
2ax4 x5 x7 x10 x217 (2x16 −
− x11 x215 x17 − x11 x15 x17 )
0
3
0
0
1
0
0
bx7
−bx6
bx9
bx8
−bx7 x11
bx7 x10
−bx6 x13
bx6 x12
bx7 x10 x17 (2x15 + 1)
bx7 x10 x16 x17 (2x15 + 1)
−bx7 x10 x15 x17 (2x15 + 1)
2bx7 x10 x217 (x16 −
− x11 x215 x17 − x11 x15 x17 )
0
0
0
Заключение. Главный результат статьи – новый алгоритм дифференцирования
функций многих переменных. Он подробно изложен в п. 3 и проиллюстрирован хоть
и модельным, но весьма содержательным и сложным по своей структуре примером
в п. 4. Теперь естественно обсудить, чем он отличается от других алгоритмов дифференцирования и в каких задачах его целесообразно применять. В различных областях
прикладной математики широко используются символьное, автоматическое и численное дифференцирование [1–3, 8–10]. Символьное дифференцирование преобразует формульное представление функции в аналогичное представление ее производной. Автоматическое дифференцирование (которое называют также дифференцированием программ) получает производную функции, тоже заданной компьютерной программой,
в виде компьютерной программы, которую далее можно применять для вычисления
значения этой функции в разных точках. Численное дифференцирование использует
приближенные формулы для производных, выведенные дифференцированием полинома, аппроксимирующего функцию, заданную в нескольких точках. Ясно, что первые
два способа позволяют, в принципе, получать значения производных достаточно точно,
и в этом их основное преимущество по сравнению с численным дифференцированием.
Главное преимущество символьного дифференцирования в том, что оно дает выражение производной в виде формулы, которую далее можно применять так, как требуется пользователю. Основной недостаток символьного дифференцирования заключается
в том, что данный метод может генерировать очень сложные выражения, содержащие много необязательных вхождений одних и тех же подвыражений. Именно потому
символьное дифференцирование эффективно в случае функций, заданных простыми
формулами, а для более сложных функций время вычисления и размер необходимой
компьютерной памяти быстро растут с увеличением сложности формульного представления этой функции.
24
Анализируя построения п. 3, несложно обнаружить, что в настоящей статье предложен новый алгоритм символьного дифференцирования, свободный от упомянутого
недостатка. Исходную функцию (систему функций) он преобразует в систему полиномов относительно исходных и дополнительно введенных переменных. При этом в результате получаются и все частные производные дополнительных переменных и исходных функций также в форме полиномов по тем же переменным. Производные более
высокого порядка исходной функции (исходной системы функций) далее можно получить дифференцированием полиномов при помощи «цепного правила».
Символьное и аналитическое дифференцирование [1–3, 9, 10] применяют, в частности, для вычисления таких величин как градиент функции, ее матрицы Якоби и Гессе
(например, в методах оптимизации). В последние годы аналитическое дифференцирование успешно используется для вычисления коэффициентов Тейлора функций, определяемых дифференциальными уравнениями [11–14]. Во всех этих задачах можно применять и предложенный в настоящей работе метод символьного дифференцирования,
а для нахождения коэффициентов Тейлора – и сам метод дополнительных переменных
(при помощи этого метода уравнения приводятся к полиномиальному виду и решаются
методом рядов Тейлора [4, 5, 15]).
Литература
1. Rall L. B. Automatic differentiation: Techniques and applications. Berlin: Springer, 1981. 171 p.
2. Naumann U. The Art of Differentiating Computer Programs: An Introduction to Algorithmic
Differentiation // SIAM. 2012. P. 340.
3. Wolfram Mathematica Documentation Center // URL: http://reference.wolfram.com/mathematica/
guide/Mathematica.html.
4. Бабаджанянц Л. К. Метод дополнительных переменных // Вестн. С.-Петерб. ун-та. Сер. 10:
Прикладная математика, информатика, процессы управления. 2010. Вып. 1. С. 3–11 .
5. Бабаджанянц Л. К., Брэгман К. М. Алгоритм метода дополнительных переменных // Вестн.
С.-Петерб. ун-та. Сер. 10: Прикладная математика, информатика, процессы управления. 2012. Вып. 2.
С. 3–12.
6. Холшевников К. В., Титов В. Б. Задача двух тел. СПб.: Изд-во С.-Петерб. ун-та, 2007. 180 с.
7. Абрамовиц М., Стиган И. Справочник по специальным функциям / пер. с англ.; под ред.
В. А. Диткина, Л. Н. Кармазиной. М.: Наука, 1979. 832 с. (Abramowitz M., Stegun I. A. Handbook of
mathematical functions.)
8. Бахвалов Н. С., Жидков Н. П., Кобельков Г. М. Численные методы. М.: Лаборатория Базовых
Знаний, 2000. 624 с.
9. Berz M., Bischof C., Corliss G.F. Computational differentiation: techniques, applications, and tools //
SIAM. 1996. 419 p.
10. Griewank A. Evaluating derivatives // SIAM. 2000. 369 p.
11. Makino K., Berz M. COSY INFINITY Version 9 // Nuclear Inst. and Methods in Physics Research.
2006. Vol. A 558, issue 1. P. 346–350.
12. Makino K., Berz M. Taylor models and other validated functional inclusion methods // Intern. J.
of Pure and Applied Mathematics. 2003. P. 239–316.
13. Abad A., Barrio R., Blesa F., Rodriguez M. Breaking the limits: the Taylor series method // Appl.
Math. and Computation. 2011. P. 7940–7954.
14. TIDES webpage URL // URL: http://gme.unizar.es/software/tides.
15. Бабаджанянц Л. К., Большаков А. И. Реализация метода рядов Тейлора для решения обыкновенных дифференциальных уравнений // Вычислительные методы и программирование. Науч.-исслед.
вычисл. центр Моск. ун-та им. М. В. Ломоносова. 2012. Т. 13. С. 497–510.
Статья рекомендована к печати проф. Е. И. Веремеем.
Статья поступила в редакцию 20 декабря 2012 г.
Документ
Категория
Без категории
Просмотров
4
Размер файла
364 Кб
Теги
алгоритм, метод, основанная, дополнительные, дифференцированный, переменных
1/--страниц
Пожаловаться на содержимое документа