ПРОГРАММИРОВАНИЕ, 2009, № 2, с. 10-27
КОМПЬЮТЕРНАЯ АЛГЕБРА
УДК 681.3.06
ПРОГРАММА НА ЯЗЫКЕ MATHEMATICA ДЛЯ "ИНТЕГРАТОРА БЕДНЯКА", РЕАЛИЗУЮЩЕГО ИДЕИ РИША И НОРМАНА*
© 2009 г. Р. Краглер
Университет прикладных наук D-88241 Вайнгартен, Германия E-mail: kragler@hs-weingarten.cle Поступила в редакцию 02.06.2008 г.
В настоящей работе путем компьютерных экспериментов исследуются достоинства и недостатки эвристического подхода, лежащего в основе метода "параллельного интегрирования". Для этой цели описывается и используется реализация данного метода на языке Mathematica, которая в ряде случаев сопоставляется также с первоначальной реализацией метода на языке Maple.
1. МЕТОД ПАРАЛЛЕЛЬНОГО ИНТЕГРИРОВАНИЯ
В своей последней (см. [1]) работе [2] Мануэль Бронштейн из INRIA (Франция) представил реализацию на языке компьютерной алгебры Maple эвристического подхода к построению неопределенных интегралов, известного в литературе как "новый метод Риша" [3], "метод Риша-Нормана" [4-6] (чаще всего) и "параллельный метод Риша" [7-9]. Представленная программа pmint названа ее автором "интегратором бедняка" (Poor Man's Integrator) и доступна в интернете [10].
В монографии [11] М. Бронштейн посвятил отдельную главу (глава 10) данному методу интегрирования, который состоит в попытке построить дифференциальное поле, содержащее как подынтегральное выражение, так и интеграл "в параллель", т.е. без последовательного (рекурсивного) построения соответствующих расширений. Фактически, метод параллельного интегрирования является не алгоритмическим, а эвристическим, поскольку он может не сработать как в теории, так и на практике при интегрировании в классе элементарных функций.
Основная идея параллельного интегрирования заключается [11] в обходе рекурсивной природы алгоритма нахождения интеграла путем рассмотрения дифференциального поля К, содержащего подынтегральное выражение / и интеграл от него, как поля рациональных функций нескольких переменных над полем (дифференциальных) констант. Тем самым, предполагается, что / 6 К, где К = ..., С = СопвЬ(К) (поле констант для К) и каждое является трансцендентным над С{Ь 1,..., 1). Как следует из теоремы Лиувилля в ее усиленной версии (см. [11], теорема 5.5.3, глава 4, стр. 145), если интеграл от подынтегрального выражения /, являющегося комбинацией элементарных функций над К, берется в этом поле (точнее, в его простом расширении), то существуют V £ К, константы с\,..., ст, принадлежащие алгебраическому замыканию поля С, и и\,..., ит 6 К, такие, что
f = Dv + Y^Cl %=1
Dm
(1)
"Перевод с англ. В.П. Гердта.
Поскольку V является отношением многочленов от ¿1,..., а щ — суть многочлены по этим переменным, метод параллельного интегрирования состоит в явном задании многочленов щ и знаменателя рациональной функции V, а также степени ее числителя. При таком задании остается произвол в коэффициентах числителя V и постоянных С{ в выражении (1).
га
U
Clear\pmint\,
pmint[f_, opt: OptionsPattern[{Extension —У 0}]] :=
Module[{ff, si, li, lin, lout, q, d, I, Iv, Id, Is, fint, lc, t, terms, list}, ff = Together[TrigToTan[f]]] list = ToIndetsAndDerivation[ff, x,t]] I ¡[list == %F ailed, Return[INT[f]]]; {ff, terms, Iv, Id} = list; q = denD[ld]] Id = q * Id;
Is = DeleteC'ases[Map[getSpecial[#, Thread[terms —У Rest[lv]]]lk, terms], Null]] ToTerms\pmIntegrate[ff, Iv, Id, q, Is, OptionValue[Extension]], terms, Rest[lv]] /. {tan —У Tan} ]
Рис. 1.
pmint := proc( f, x)
local ff, si, li, lin, lout, Id, q, d, I, vars, dx, Is, fint, lc] ff := eval(convert(f, tan))] ] ф convert trigs to tan
si := select(proc(d) diff(d,x) <> 0 end, indets(f/)); si := select(proc(d) diff(d, x) <> 0 end, indets(map(diff, si, x))) union si;
H := [op(si)]; ф list of terms in integrand and its derivative
lin := [seq(d = Hools/genglobal1 (x), d = lin)]; ф substitution terms — > indets
lout := [seq(rhs(d) = lhs(d), d = lin)]; ф substitution indets — > terms
Id := subs(lin, map(diff, li, x))] ф derivatives of all the terms
q := lcm(seq(denom(d), d = Id))] ф denominator of the total derivation
I := [seq(normal(q * d), d = Id)]] ф q * derivatives of all the terms
vars := map(lhs,lout)]
dx := totalDerivation(vars, l)] ф vector field dx = q * d/dx
Is := [seq(getSpecial(d, lin), d = li)]] ф list of known Darboux for dx
fint := subs(lout, pmlntegrate(subs(lin, f f), dx, q, vars, Is))]
lc := select(proc(d) converted, string)[l] =" _"; end, indets(fint, name))]
subs({seq(d = 0, d = lc minus si)}, fint)]
end]
Рис. 2.
Тем самым, задача сводится к решению системы линейных уравнений, вытекающих из равенства (1), относительно постоянных сг- и коэффициентов числителя рациональной функции V. Решение получаемой линейной системы осуществляется стандартными методами линейной алгебры.
Если система линейных уравнений имеет решение, то интеграл от / находится сразу же очевидным интегрированием выражения (1). В то же время, отсутствие (несуществование) реше-
ния, вообще говоря, не означает, что интеграл от / не берется в рассматриваемом классе функций. Отсутствие решения может просто означать, что задание указанной структуры является неверным. Все известные в литературе подобные "параллельные" подходы к интегрированию обладают таким недостатком, и для каждого из них имеются примеры подынтегральных выражений в классе элементарных функций, которые допускают интегрирование в этом классе, но интегралы от которых не находят-
Clear\pmlntegrate];
pmlntegrate[f _,lv_List,ld_List,q_,ls_: {},ext_: 0]:=
Module[{splq, s, df, spl, cden, dg, x, monomials, cand, lunk, sol, i, AA}, DRPrint["Visit pmlntegrate."]; splq = splitF actor[q,lv ,ld\, s = First[splq\,
Scan[If[Last[#],s = s * First[#\]k, /s];
x = First\lv\,
df = Denominator[f\,
spl = splitFactor[df,lv,ld];
cden = s * First[spl] * deflation[Last[spl], Iv, ld\, dg = 1 + ioia/Deg[s, Iv] +
Ma«[ioia/_De(jf[#Mmeraior[/], /i>], totalDeg[Denominator[f], Iv]]; monomials = enumerateMonoms[lv, dg]; DRPrint["There are", Length[lv],
"new variables in the integrand and the guess bound for deg is", dg, "therefore the number of monomials is", Length[monomials]," ."];
If[Length[monomials] > 800, DRPrint["\n Since the number of monomials is bigger than 800, then
'pmint' failed. It is possible that the integral is doable,"]; DRPrint["but it will require too much time. The bound 800 can be changed to a smaller one"]; Return[INT[f]]]; cand = Total[Table[AA[i] * monomials[[i]],
{i, Length[monomials]}]]/cden; lunk = Table[AA[i], {i, Length[monomials]}]; sol = trylntegral[f, Iv, Id, q, cand, lunk,
First[spl], (* normal part of df *)
Last[spl], (* special part of df *)
Last[splq], (* special part of q *)
/ S у & tie i j ^
If[First[sol], INT[f], Last[sol]] ]
Рис. 3.
ся "параллельным" методом. Это свидетельствует об эвристическом характере данного метода. Тем не менее, метод допускает сравнительно (по отношению к полному алгоритму Ри-ша) легкую программную реализацию на языках систем компьютерной алгебры, таких, как Maple и Mathematica, и довольно хорошо работает на практике. Ниже мы продемонстрируем эти преимущества рассматриваемого метода на конкретных примерах.
Реализация на языке Maple эвристического параллельного интегрирования в виде программы pmint была выложена М. Бронштейном в Интернет в мае 2005 года [10] и содержит всего лишь 95 строк кода. Данная программа ориен-
2. ПРОГРАММА НА ЯЗЫКЕ MATHEMATICA ДЛЯ ПАРАЛЛЕЛЬНОГО ИНТЕГРИРОВАНИЯ
pmlntegrate := proc(f, d, q, vars)
local Is, splq, s, ff, df, spl, cden, dg, monomials, cand, lunk, sol, i; if nargs = 5 then Is := args[b]; else Is := []; /г; splq := split F act or (q, d); s := splq[ 1];
for i to nops(ls) do if /s[i][2] then s := s * /s[i][l]; /г; od; ff := normal (f); df := denom(ff); spl := splitFactor(df, d); cden := s * sp/[l] * de//aiion(sp/[2], d);
dg := 1 + degree(s) + max(degree(numer(ff)), degree(denom(f /)));
monomials := [op(enumerateMonoms(vars, dg))];
cand := add('_A'[i] * monomia/s[i], г = 1 ..nops(monomials))/cden]
lunk := {seg('_=4'[i], г = 1 ..nops(monomials))}]
sol := trylntegral(f, d, q, vars, cand, lunk,
spl[l],spl[2],splq[l],ls,0y, if sol[l] then sol := trylntegral(f, d, q, vars, cand, lunk,
spl[l],spl[2],splq[l],ls, I); fi; if sol[l] then Int(f); else sol[2]; fi; end;
Рис. 4.
тирована на интегрирование трансцендентных функций, как элементарных, так и специальных.
Л. Медина и А. Павлык разработали [12] концепцию prriint для системы Mathematica. Во введении мы подчеркивали, что задача "перепрограммирования" prriint с языка Maple на язык Mathematica является нетривиальной из-за различной природы и функциональности встроенных системных функций. Автором настоящей работы произведена доработка программы prriint и выполнено ее интенсивное тестирование на различных примерах для версии 6.0.2 системы Mathematica.
Ниже мы представим код prriint на языке системы Mathematica в виде отдельных процедур и, где это возможно, сравним его с кодом соответствующих процедур в Maple.
(1) Процедуры: pmint, pmlntegrate, trylntegral, getSpecial, rriyFactors, enurrierateM ononis, splitFactor, deflation.
Процедура pmint на языке Mathematica представлена на рис. 1.
Для сравнения приведем соответствующую процедуру pmint на языке Maple (см. рис. 2).
На рис. 3 приведена процедура pmlntegrate на языке Mathematica.
Соответствующая процедура pmlntegrate на языке Maple показана на рис. 4.
Процедура trylntegral на языке Mathematica представлена на рис. 5.
Соответствующая процедура trylntegral на языке Maple изображена на рис. 6.
Процедура getSpecial на языке Mathematica (для получения известных полиномов Дарбу) показана на рис. 7.
Соответствующая процедура getSpecial на языке Maple возвращает многочлены Дарбу (см. рис. 8).
Процедура myFactors на языке Mathematica приведена на рис. 9.
Соответствующая процедура myFactors на языке Maple представлена на рис. 10.
Процедура enumerateMonoms на языке Mathematica вычисляет все мономы от переменных "lv" такие, что их полная степень не больше "deg" (см. рис. 11).
На рис. 12 изображена соответству
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.