Алгоритм дифференцирования, основанный на методе дополнительных переменных Текст научной статьи по специальности «Математика»

УДК 517.9:519.6 Вестник СПбГУ. Сер. 10, 2013, вып. 2

К. М. Брэгман

АЛГОРИТМ ДИФФЕРЕНЦИРОВАНИЯ, ОСНОВАННЫЙ НА МЕТОДЕ ДОПОЛНИТЕЛЬНЫХ ПЕРЕМЕННЫХ

Введение. Средства символьного и автоматического дифференцирования функций многих переменных [1,2] настолько важны и востребованы, что входят во все основные математические пакеты компьютерной алгебры (см., например, [3]). В настоящей статье предлагается новый алгоритм нахождения производных системы функций многих переменных весьма общего вида. Речь идет о функциях, которые можно получить из своих аргументов при помощи конечного числа операций +, —, х,/ и суперпозиций функций, удовлетворяющих полным полиномиальным системам уравнений в частных производных (и, в частности, системам обыкновенных дифференциальных уравнений - ОДУ). Класс таких функций включает в себя элементарные функции, специальные функции математической физики (за исключением теоретико-числовых функций) и многие другие. Алгоритм основан на методе дополнительных переменных (МДП) и использует библиотеку, содержащую набор таких функций (а точнее, полиномиальных дифференциальных уравнений для этих функций). Как и в статьях [4, 5], предполагается, что исходная система функций и функции библиотеки могут зависеть также и от конечного числа параметров. Последовательное введение новых (дополнительных) переменных сводит все функции исходной системы к полиномиальным функциям по исходным и новым аргументам, причем новые аргументы удовлетворяют полной полиномиальной системе уравнений в частных производных (по исходным аргументам). Получение производных исходной системы сводится после этого к использованию «цепного правила» дифференцирования сложных функций и, в итоге, к перемножению и сложению полиномов, а результаты выражаются полиномами по исходным и новым переменным. Для удобства читателя в начале статьи (в п. 1) сразу приводится нетривиальный пример, поясняющий предлагаемый алгоритм и его отличие от известных алгоритмов символьного и автоматического дифференцирования: рассматриваются две функции трех переменных и пяти параметров и показывается, как после введения нескольких дополнительных переменных вычисление частных производных любого порядка для этих функций сводится к сложению и умножению полиномов. В п. 2 приводятся необходимые сведения о классах функций, дифференциальных уравнениях и библиотеках функций, в п. 3 описывается сам алгоритм, а в п. 4 - пример его применения.

Отметим наконец, что для нахождения производных функций многих переменных мы применяем те же средства, что и в статье [5], поэтому здесь, насколько возможно, используем аналогичные обозначения, а также аналогичные последовательность и способ изложения. Это позволило предоставить пользователю возможности преобразования дифференциальных уравнений к полиномиальной форме и дифференцирования функций в рамках схожих алгоритмов метода дополнительных переменных. Естественно, что читателю рекомендуется сначала ознакомиться со статьей [5] (хотя это и не обязательно), а затем уже с настоящей работой, содержащей алгоритм еще одного применения МДП.

Брэгман Константин Михайлович — аспирант, Санкт-Петербургский государственный университет; e-mail: [email protected].

© К. М. Брэгман, 2013

1. Поясняющий пример. Рассмотрим систему из двух функций трех переменных уод = х-1 (ах\ + Ьх2) хз 81П(^ж2),

Уо,2 = 1п12 Х1(с3Х1Х2х3 совз(их2) + соШ(с^) 81п(с^х2))3, где х1,х2, хз - переменные, а а, Ь, с,й,ш - параметры. Вводя еще переменные

• / \ / \ 1 —1 Х3 СОв(шХ2)

<1 = Я1п(^х2), <2 = СОВ(^х2), <з = 1п х1, <4 = х1 , <5 = х1 ,

вместо уо,1, уо,2 получаем полиномы (по переменным х = (х1, х2, хз), < = (<1,..., <£>5))

У1 = (ах1 + Ьх2)7хз<1<5, у2 = <32(с3х1х2хз<3 + соШ(с^ ¿)<1)3,

причем < является решением полной системы:

д<1/дх1 = 0, д<1 /дх2 = ш<2, д<1/дхз = 0, д<2/дх1 = 0,

д<2/дх2 = -ш<1, д<2/дхз = 0, д<з/дх1 = <4, д<з/дх2 = 0, д<з/дхз = 0, д<4/дх1 = -<\, д<4/дх2 = 0, д<4/дхз = 0, д<5/дх1 = хз<2<4<5, д<5/дх2 = -ихз<1<з<5, д<5/дхз = <2<3<5. Если использовать обозначения

^уо,г = д<1+<2+<3уо,г/дх1 дх%дхЗ3, г = (¿1,12,13), Н € [0 :

е1 = (1, 0, 0), е2 = (0, 1, 0), ез = (0, 0, 1),

то для вычисления искомых производных можно использовать рекуррентные формулы («цепное правило» дифференцирования сложных функций)

Оек уо,г = дуг/дхк дуг/д<т ■ д<т/дхк

1

&+екуо,г = дБ'уо,г/дхк + ^ дБ'уо,г/д<т ■ д<т/дхк, г =1, 2; к = 1, 2, 3,

т=1

причем каждая из получаемых величин Векуог,...,Вг+екуог будет полиномом по х1,х2, хз, <1,...,<5, коэффициенты которого зависят от а = (а, Ь, с, ¿, и).

Для вычисления производной дг1+г2+®3уо ,г/дх1 дх2,2дх^З для любых значений х,а можно сначала определить <, а затем значение полинома В1 у о , г, причем его можно рассчитать по явной формуле (если она получена заранее по приведенным выше рекуррентным соотношениям) или непосредственно по рекуррентным соотношениям.

Таким образом, можно сказать, что предлагаемый алгоритм использует, как и общепринятые методы символьного и автоматического дифференцирования, «цепное правило», но благодаря технике МДП (величины дуг/дхк, дуг/д<т, д<т/дхк в этом случае оказываются полиномами по исходным и дополнительным переменным) алгоритм применения этого правила сводится сначала к вычислению дополнительных переменных и затем к перемножению и сложению полиномов.

Заметим, что алгоритм не отвечает на два ключевых вопроса: как автоматизировать выбор дополнительных переменных и как найти для них полную систему уравнений в частных производных с полиномиальными правыми частями. Ответы на эти вопросы будут даны в п. 3.

2. Дифференциальные уравнения, классы функций, библиотеки. Здесь приводятся необходимые сведения из статей [4, 5]. В п. 2.1 вводятся в рассмотрение полные полиномиальные системы дифференциальных уравнений и классы функций, которые таким системам удовлетворяют; в п. 2.2 - классы функций многих переменных, на которые ориентирован предлагаемый в настоящей статье алгоритм (см. п. 3). В п. 2.3 обсуждается основной используемый в этом алгоритме инструмент - библиотека.

2.1. Дифференциальные уравнения, классы £а. Положим

х = (х\,...,хт) е Ят, £ = (¿1,...,^) € Я8, а =(а1,...,аш) € Яш

и рассмотрим систему уравнений

дхг/дг^ = Ц(х,а), I е [1 : т], ] е [1 : в],

где все функции - полиномы по х1,..., хт с коэффициентами, зависящими от а. Такую систему называют полиномиальной. В случае в = 1 это система ОДУ.

Говорят, что скалярная функция р аргумента х = (х1,...,ха) удовлетворяет полиномиальной системе, если она является одной из компонент вектор-функции (того же аргумента х = (х1,... ,ха)) - решения некоторой полиномиальной системы. Класс скалярных функций аргумента х = (х1,.. .,ха), удовлетворяющих полиномиальной системе, обозначим £а. Ясно, что £1 С £2 С ... . Все элементарные функции и большое число специальных функций принадлежат £1 (и тем более £а при а > 1).

2.2. Классы Тт. Классом Тт, а, т е [1 : называют множество скалярных функций аргумента х = (х1,... ,хт), которые можно получить из х1,... ,хт при помощи конечного числа операций +, —, х,/ и конечного числа функций р1,...,р; е £а и их суперпозиций.

2.3. Расширения и библиотеки. Если р е £а, то существует вектор-функция (р1, . ..,рп) - решение полиномиальной системы, где р1 = р. Расширением р называют множество {р1, ...,рп}. Объединение расширений нескольких функций называют библиотекой. Фактически, библиотека содержит набор полиномиальных систем, а не задач Коши. Всякое подмножество библиотеки называют подбиблиотекой, если оно само является библиотекой. Библиотеку называют автономной или неавтономной соответственно тому, какими уравнениями она представлена - автономными или неавтономными. Естественно пользоваться как автономными, так и неавтономными библиотеками, хотя можно было бы ограничиться автономными библиотеками, так как неавтономную библиотеку всегда можно свести к автономной.

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

имевшуюся нумерацию разделов можно сохранить, а новым разделам (если появятся) присвоить новые номера (старые разделы могут пополниться новыми функциями).

Для пользователя библиотека - пополняемая таблица, состоящая из потенциально неограниченного количества строк и одиннадцати столбцов, полное описание которых дано в п. 2.1 статьи [5].

3. Алгоритм дифференцирования. Рассмотрим систему функций

yr = fr(x), r е [1: N] (или y = f (x)), (1)

при x = (xi,...,xm) е Rm, y = (yi,...,yN), f = (fi,...,fN) е RN и при условии, что все fr (x) принадлежат классу F^ и записаны в терминах функций имеющейся в нашем распоряжении библиотеки. Кроме того, будем считать, что fr (x) могут зависеть от параметров а = (ai, ...,аш) е Rw, хотя мы и не указываем на это явно в обозначениях. Мы собираемся здесь вывести формулы для производных dl1+ +imyr/dxl ... dx'm.

В п. 3.1 покажем, как, вводя дополнительные переменные при помощи библиотеки, преобразовать систему к полиномиальной форме, и выведем рекуррентные формулы для производных введенных переменных по xi,..., xm. В п. 3.2 получим формулы для искомых производных, а в п. 3.3 опишем схему алгоритма.

3.1. Элементарное преобразование системы. Преобразование исходной системы к полиномиальной состоит из конечного числа K последовательных шагов, причем на каждом из них проводится однотипное «элементарное» преобразование, которое мы далее рассмотрим. Ради симметрии обозначений, для величин yr,y,fr,f,x,m в (1) будем использовать также и обозначения yo,r,yo,fo,r,fo,x0,m(0) соответственно. На каждом шаге вводятся дополнительные переменные, причем в новых переменных система (1) представляется в виде (к - номер шага)

yk,r = fk,r(xk), r е [1: N] или y = fk(xk), к е [1 : K], (2)

x (x ; X ) (x1, ..., xm(k-1); xm(k-1) + 1,..., xm(k- 1)+n(k) ) (x1,..., xm(k) ),

m(k) = m(k - 1) + n(k), к е [1 : K], x0 = x = (xb...,xm(o)), m(0) = m, yo,r = Уг , fo,r = fr,

а натуральные числа n(k) и функции f1jr, f2,r,... определяются на каждом очередном шаге. При к = 0 системы (1) и (2) совпадают. Первый шаг состоит в переходе от системы (1) к системе (2) при к = 1, а на последнем шаге, при к = K, система (2) становится полиномиальной. Перейдем к описанию к-го шага (к е [1 : K]).

а). Ищем в выражениях fk-1,r(xk-1), r е [1 : N] функцию вида ф^Р k(xk-1))p(k), где Рk (xk-1) = (PM(xk-1 ),...,Pk,#(k) (xk-1)), Pk,v - алгебраические полиномы, фk = ф^р-i^,... ,p$(k)) - функция из библиотеки, а р(к),д(к) е [1 : +го).

б). Пусть фk = у>1, а ^1,..., ^n(k) - все функции расширения функции ф k и пусть

d^i/dpv = QvKl(pk; fk), l е [1: п(к)], ^ е [1 : д(к)], (3)

где pk = (p1,.. .,p$(k)), <fk = (<f1,.. .,^n(k)), а Qvkii(pk; ¥k) - полиномы по всем п(к) + $(к) своим аргументам. Вводим п(к) новых дополнительных переменных

xm(k-1)+1 = ^1(Рk (xk-1)),..., xm(k) = Pn(k) (Pk (xk-1)). (4)

в). Заменяя в /к-1,г (хк г е [1 : М], все функции р1(Рк (хк 1 )),...,рп(к)(Р к(хк х)) на новые переменные хт(к-1)+ъ ..., хт(к) соответственно, получаем

Ук,т = 1к,т(хк), г е [1: М]. (5)

г). При I е [1 : и(к)\, в е [1 : т], выписываем производные (см. (3), (4))

причем ¿/¿х8 - производная сложной функции по аргументу х8, в е [1 : т] при постоянных значениях аргументов х^, ] е [1 : т], ] = в, а хп(хк-1) = хп при п е [1 : т].

Для того чтобы воспользоваться формулами (6), (7) для последовательного (при к = 1,...,К) нахождения производных дополнительных переменных по аргументам х1,..., хт, необходимо на каждом шаге иметь в своем распоряжении величины

n(k), Qk,i(pi>--->p#(h); ), pk,v (xk 1),

l e [1 : n(k)], V e [1 : tf(k)], '

которые уже известны (см. a), б)).

3.2. Производные исходной системы. Если использовать обозначения

DV = fr(x)/dx1 ...dx%, i =(ii,...,im), ij e [0:+rc),

ei = (1, 0,...,0),...,em = (0,...,0,1), то, применяя «цепное правило» дифференцирования, получаем

д^о д^п dLX о

j M=m(0)+1 М j

ni+e, дРгуг дРгуг dxм

D Уг = — + -ЦГ'^ rell.Nl 3G[l.m],

j M=m(0)+1 M j

причем производные dx^/dxj вычисляются по формулам (6), (7) и dxм/dxj = S^j при e [1 : m].

Найдем более удобный вариант этих формул, для чего полиномы yr в (5) при k = K запишем в виде

u

yr Ху ar,nx П , x (x1, ..., xm (K)),

n=0 (9)

Kn) = (\,1,..., (k) ), Г e [1 : N],

где K(0) = (0,...,0), xA(1) = x1,...,xA(m(K)) = xm(K), а xA(m(K)+1),...,xx(u) - все различные нелинейные мономы, из которых составлены полиномы fxr, r e [1 : N].

Полагая В(0'"''0)уг = уг, обобщим представление (9) на В1 уг:

«(«)

ВУ = ^ аг,п (г)хХ(г'п), X = (хЬ ...,Хт(К)),

п=0

Л(г, ^ = (Лп,1(«),..., \,т (К)(^)), г е [1 : N].

Обозначения здесь аналогичны и очевидны, причем эта формула включает в себя (9) как частный случай при г = (0,..., 0). Дифференцируя ее, имеем следующие рекуррентные соотношения для нахождения искомых производных:

u(i) t m(K)

d>^c / i

D^'Vr = E • -¿Г I • (10)

n=o \ /=1

Вычисления по формулам (10) могут отличаться тем, в какой последовательности находятся производные (полиномы по исходным и введенным аргументам). Производные dxм/dxj находятся заранее по рекуррентным формулам (6), (7).

3.3. Схема алгоритма.

Шаг 1. Заносим систему (1) в таблицу (например, в пакете Wolfram Mathematica). Введенная система может содержать параметры. Если для обозначения каких-то из этих параметров использованы символы вида aY, где y - натуральное число, то полагаем, что Q равно максимальному из таких y, а если таких параметров нет, то полагаем Q = 0. Параметры, от которых зависят функции подбиблиотеки, переобозначим символами i, 2 ,••• и запомним таблицу соответствий новых и старых обозначений. Образуем подбиблиотеку, состоящую из объединения расширений всех функций, участвующих в написании исходной системы.

Шаг 2. Преобразуем систему шаг за шагом по алгоритму п. 3.1, получая на каждом из них систему вида (5), пока она не окажется полиномиальной.

Шаг 3. Вычисляем производные по формулам (6), (7).

Заметим, что в программе можно эти производные вычислять и последовательно в рамках шага 2.

Шаг 4. Рассчитываем производные исходной системы по алгоритму п. 3.2.

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

4. Пример. Чтобы ограничиться той же библиотекой, что и в статье [5], рассмотрим систему из шести функций, которые немногим отличаются от правых частей полной системы уравнений, использованной там в качестве примера:

yi = fi(xi, Х2, хз) = x3(sincos(aln2 X2 + 6x3) + bln4 X2) + + sin5(aln2 x2 + bx3) + Dv(g2,1, —1; x1), У2 = f2(xi,x2,x3) = (sin cos(a ln2 x2 + 6x3) +b ln4 x2) + + sin5(aln2 x2 + bx3) + Dv(g2,1, —1; xi), Уз = f3(xi,x2,x3) = x2 cossin(aln x2 + bx3) + + cos4(aln2 x2 + bx3) + EK2(sinsin(aln2 x2 + bx3), xi),

У4 = f4(xi, x,2, x3) = x2 cossin(a ln2 x2 + bx3) + + cos4(aln2 x2 + bx3) + EK(sinsin(aln2 x2 + bx3), xi),

(11)

Вход и заведение системы Сообщение об ошибке

Нет

Возможно ли построение подбиблиотеки

Да

Построение подбиблиотеки для исходной системы. Формирование разделов и таблицы параметров

Заведение и проверка данных

Подготовка

Является ли полиномиальной полученная система?

Нет

Найти в системе

функцию, аргументы

которой - полиномы

Преобразо-

вания

Вывод данных

Блок-схема алгоритма дифференцирования

У5 = /б(жьЖ2,жз) = x2(ch2(aln2 X2 + 6x3) + + sin(a ln2 x2 + bx3)) + sh5(a ln2 x2 + bx3), У6 = /e(xi,x2,x3) = xi(ch (aln x2 + 6x3) + + sin(aln2 x2 + bx3)) + sh5(aln2 x2 + bx3),

где a = (ai, a2, аз) = (a, b, g) - вещественные постоянные (параметры).

Упомянутая выше библиотека является объединением шести подбиблиотек (Пб), которые мы сейчас рассмотрим вместе с соответствующими дифференциальными уравнениями. Функции обозначаются p,pi,..., а их аргументы - p,pi,... [5].

Пб1. Состоит из одной функции одного аргумента - p = inv(p) = p-i, которая удовлетворяет уравнению dp/dp = — p2, а расширение этой функции состоит из нее одной.

Пб2. Состоит из двух функций одного аргумента - pi = lnp, р2 = inv(p), которые удовлетворяют системе dpi/dp = р2, dp2/dp = —р2, причем расширение функции pi состоит из функций pi,p2, а расширение функции р2 - из нее одной (т. е. совпадает с Пб1).

Пб3. Состоит из двух функций одного аргумента - pi = sinp, p2 = cosp, которые удовлетворяют системе dpi/dp = p2, dp2/dp = —pi, а расширение каждой из функций pi, p2 состоит из функций pi,p2.

Пб4. Состоит из двух функций одного аргумента - pi = sh(p), p2 = ch(p), которые удовлетворяют системе dpi/dp = p2, dp2/dp = pi, а расширение каждой из функций pi, p2 состоит из функций pi,p2.

Пб5. Состоит из четырех функций двух аргументов - pi = EK(pi,p2) = E(e,M), где EK = E - эксцентрическая аномалия, pi = e - эксцентриситет, p2 = M - средняя аномалия (см. [6]),

p2 = EKs(pi,p2) = sin pi, p3 = EKc(pi,p2) = cos pi, p4 = EKi(pi,p2) = (1 — pi cospi)-i, которые удовлетворяют системе

dpi/dpi = p2p4, dpi/dp2 = p4, dp2/dpi = p2p3p4, dp2/dp2 = p3 p4,

dp3/dpi = —p2p4, dp3/dp2 = —p2p4,

dp4/dpi = p3 p2 — pip2p4, dp4/dp2 = — pip3p2,

причем расширение функции pi состоит из pi,p2,p3,p4, а расширение каждой из функций p2,p3,p4 - из всех этих трех функций.

Пб6. Состоит из двух функций pi, p2 аргумента p, причем первая из них - функция Вебера Dv, а вторая - ее производная, которая обозначена DV. Как известно [7], они удовлетворяют системе уравнений

dpi/dp = p2, dp2/dp = — (ap2 + bp + c)pi,

где a,b, c - параметры. Расширение функций pi, p2 состоит из функций pi ,p2.

Переходим к реализации алгоритма, описанного в п. 3.

Шаг 1. Пользователь заводит систему (11) (например, в пакете «МаШетайса»).

Шаг 2. Преобразуем систему шаг за шагом по алгоритму п. 3.1 (шаг 2.1, шаг 2.2,...), получая на каждом из них систему вида (5), пока она не окажется полиномиальной. Всего шагов оказалось семь (т. е. К = 7). Ради экономии места, выпишем сразу все дополнительные переменные, расположив их по строкам так, что первая из переменных в к-й строке соответствует функции фк = Фк(р1, ■ ■ ■ ,Р${к)), найденной на к-м шаге, а все переменные в этой строке соответствуют функциям ее расширения фк = уд,..., уп(к) (см. п. 3.1 а), б)):

ж 14 = EK(xii, xi)

Х4 = ln Ж2, Х6 = sin(ax4 + Ьжз), Х8 = sh(ax4 + bxs), xio = cos x6, xi2 = sin x7, , xi5 = EKs(xii,xi), xis = Dv(g2,1, -1; xi),

x5 — x2 ,

x7 = cos(ax4 + bx3), xg = ch(ax4 + bx3),

xii = sin x6, (12)

xi3 = cos x7,

xi6 = EKc(xii,xi), xi7 = EKi(xii ,xi), xig = DV(g2,1, -1; xi).

После седьмого шага исходная система (11) приобретает полиномиальную форму:

V7,i = Í7,i(xi, . .., xig ) -У7,2 = /7,2 (xi, . . . ,xig)

x3 (xi2 + bx4) + x5 + xis, = (xi2 + bx4) + x5 + xis,

У7,3 = /7,3 (xi, . . . , xig ) = x2xio + x7 + xM, У7,4 = /7,4 (xi, . . . ,xig) = x2xio + x7 + xi4, У7,5 = /7,5 (xi, . . . , xig ) = x2(xg + xe) + xS, У7,6 = /7,6(xi, . . . ,xig) = xi(xg + xe) + xl|.

Шаги 3, 4- Для того чтобы воспользоваться формулами (6), (7) для последовательного (при к = 1,..., 7) нахождения производных x4,..., xig по xi, x2, x3, необходимо на каждом шаге иметь в своем распоряжении величины (8). Пользуясь библиотекой, при пошаговом введении переменных (а они у нас так и записаны в (12)), их можно получить также по шагам. Последовательность действий опишем на примере шестого шага, на котором найдена функция ф6(P6д(x5),P6,2(x5))2 = (EK(xii, xjJ)2.

Как видим, P6,i(x5) = xii, P6,2 (x5) = xi, а ф6(рьр2) = EK(pi,p2) - функция из подбиблиотеки Пб5. Так как расширение ф6(рьр2) = y1(p1,p2), y2(p1,p2) = sin yi, у3 (pi, p2) = cos у i, ^4 (p i ,p2) = (1 - pi cos у i )_1 функции ф6 (pi, p2) состоит из четырех функций двух аргументов, то n(6) = 4, $(6) = 2, а так как функции у i ,...,у4 удовлетворяют системе (Ql6J = Q^ (pip у i, ...,у 4))

dipi dpi

С!.. ~~ ^бД

У2У4,

chp1 dp2

Cl.. ~~ Ql,í

у 4,

ду 2 dpi

о.. ~~ Qe,2

У2У3 У4,

ду 2 dp2

a.. ~~ Qe,2

У3У4,

д(рз dpi

i2 = Q6,3 =

дñ = Q6'3 =

д(Р4 _ ni дР1 ~ Уб'4

у3у24

23

P1У2У3,

д(Р4 _ П2 др2 Q6A

-P1У4У2,

то и полиномы Qvk 1 известны. Величины (8) для всех шагов получаются аналогично, они представлены в табл. 1. В табл. 2 приведены используемые в формуле (7) частные производные дРк,и (хк-1 )/дхп, и, наконец, в табл. 3 даны рассчитанные по формулам (6), (7) выражения (полиномы) для ¿х^/¿х^, ¡1 = 1,..., 19; ] = 1, 2, 3.

Теперь по формулам (10) можно получить выражения для производных функций системы (11) по Ж1,Ж2,хз. Пусть, например, требуется многократно вычислять производные дув/Эх^, ] = 1, 2, 3. Выводя формулы для этих величин, имеем

дуб/дх1 = (х2 + Хб) + Х1 (2хд • ¿х^/ёх! + 1хв/1х 1) + 5x8 ' 1хв/1х1 = х?д + хв,

дуб/дх2 = 1х1/1х2 • (х9 + хв) + х1 (2хд • 1хд/1х2 + 1хв/1х2) + 5х^ • 1х8/1х2 = = 4ах1х4х5х8хд + 2ах1х4х5х7 + 10ах4х5 х3хд = 2ах4х5(2х1х8 хд + х1х7 + 5х§хд), дуб/дхз = 1х1/1хз • (хд + хв) + х1 (2хд • 1хд/1хз + 1хв/1хз) + 5х4 • 1х8/1хз =

= х1(2Ьх8хд + Ьхч) + 5Ьх^д.

Таблица 1. Величины $(к), п(к), Рк^Рь •••> Р^(к); ,---,Рп(к)), Рк,^(хк_1)

к ■в (к) п(к)

1 1 2 Р1,1 = х2 Я\Л=<Р2, <9\,2 = ~Ч>2

2 1 2 Р2 д = ах\ + Ьж3 $2,1 = ¥>2, <52,2 = -VI

3 1 2 Рзд = ах\ + Ьж3 <9зд = Ч>2, <9з,2 = VI

4 1 2 -Р4Д = хв <94,1 =-^2, <94,2 =

5 1 2 -РбД = <9БД=¥>2, <9Б,2 = -VI

6 2 4 Лзд = XII, Рб,2 = XI <9бД = ¥>2¥>4, <9б,1 = <9 6 2 = ¥>2¥>з¥>4, <9 6,2 = ¥>3¥>4, <9б,з = -¥>¡¥>4, <926,3 = -^2^4, <9б,4 = ¥>3¥>1 <91,4 =

7 1 2 Рт,1 = XI <97,1 = ¥>2, <97,2 = С1 - 92Р2 ~Р)<Р 1

Таблица 2. Ненулевые величины 9Рк^ (хк 1)/дхп

к 6> Рк.Дхк-1)/^,,

1 №1,1/^2 = 1

2 ЭР2,1 /дхз = Ь, дР2,1 /дх4 = 2ах4

3 дРз х/дхз = Ь, дРз х/дх4, = 2ах4

4 дР41х /Эхе = 1

5 дРб,1/дх7 = 1

6 9Рб,1/дхц = 1, дРе^/дх! = 1

7 дР7Л/дх-1 = 1

Таблица 3. Производные dxм/dx ., ^ = 1,..., 19; ] = 1, 2, 3

м 3

1 2 3

1 1 0 0

2 0 1 0

3 0 0 1

4 0 Хб 0

5 0 -х\ 0

6 0 2ах4хьх7 Ьх7

7 0 — 2ах4хвхв -Ьх6

8 0 2ах4хвхд Ьх 9

9 0 2аХ4ХдХ8 Ьх8

10 0 — 2ах4хьхгхц —Ьх7хц

11 0 2ах4хцх7хю Ьх7х ю

12 0 — 2а,Х4ХьХбХ1з —Ьх 6X13

13 0 211X4X5X6X12 Ьх&х 12

14 а; 17 2ах4х5х7хюх17(2х15 + 1) ЬЖ7Ж10Ж17(2Ж1Б + 1)

15 Ж161С17 2ах4хвх7х10х1вх17(2х1в + 1) Ьх7х10х1вх17(2х1В + 1)

16 -Ж1БЖ17 -2аж4ЖБЖ7ЖюЖ1БЖ17(2ж1Б + 1) -ЬЖ7Ж10Ж1БЖ17(2Ж1Б + 1)

17 -хцх1бх^ 2ах4х5хухю%1у(2х1в — - ХЦХ^Х1Т - ЖЦЖ1БЖ17) 2Ьх7х10х(7(х16 - - ХЦХ^Х 17 - ЖЦЖ1БЖ17)

18 Х19 0 0

19 (1 -д2х\ - ~ Ж1)Ж18 0 0

Заключение. Главный результат статьи - новый алгоритм дифференцирования функций многих переменных. Он подробно изложен в п. 3 и проиллюстрирован хоть и модельным, но весьма содержательным и сложным по своей структуре примером в п. 4. Теперь естественно обсудить, чем он отличается от других алгоритмов дифференцирования и в каких задачах его целесообразно применять. В различных областях прикладной математики широко используются символьное, автоматическое и численное дифференцирование [1-3, 8-10]. Символьное дифференцирование преобразует формульное представление функции в аналогичное представление ее производной. Автоматическое дифференцирование (которое называют также дифференцированием программ) получает производную функции, тоже заданной компьютерной программой, в виде компьютерной программы, которую далее можно применять для вычисления значения этой функции в разных точках. Численное дифференцирование использует приближенные формулы для производных, выведенные дифференцированием полинома, аппроксимирующего функцию, заданную в нескольких точках. Ясно, что первые два способа позволяют, в принципе, получать значения производных достаточно точно, и в этом их основное преимущество по сравнению с численным дифференцированием. Главное преимущество символьного дифференцирования в том, что оно дает выражение производной в виде формулы, которую далее можно применять так, как требуется пользователю. Основной недостаток символьного дифференцирования заключается в том, что данный метод может генерировать очень сложные выражения, содержащие много необязательных вхождений одних и тех же подвыражений. Именно потому символьное дифференцирование эффективно в случае функций, заданных простыми формулами, а для более сложных функций время вычисления и размер необходимой компьютерной памяти быстро растут с увеличением сложности формульного представления этой функции.

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

Символьное и аналитическое дифференцирование [1-3, 9, 10] применяют, в частности, для вычисления таких величин как градиент функции, ее матрицы Якоби и Гессе (например, в методах оптимизации). В последние годы аналитическое дифференцирование успешно используется для вычисления коэффициентов Тейлора функций, определяемых дифференциальными уравнениями [11-14]. Во всех этих задачах можно применять и предложенный в настоящей работе метод символьного дифференцирования, а для нахождения коэффициентов Тейлора - и сам метод дополнительных переменных (при помощи этого метода уравнения приводятся к полиномиальному виду и решаются методом рядов Тейлора [4, 5, 15]).

Литература

1. Rail 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 г.