Mathematica Implementation of M



Programming and Computer Software (2009) Vo. 35, No. 2, pp 63-78

Mathematica Implementation of M. Bronstein's Poor Man's Integrator

Parallel Integration Method – Risch Norman Algorithm

Prof. Dr. Robert Kragler, Weingarten Univ. of Applied Sciences, Germany

email : kragler@hs-weingarten.de

web :

Abstract :

In May 2005 Manuel Bronstein presented on the same CAS seminar in Dubna his implementation of the so-called "Poor Man's Integrator". The method of "Parallel Integration" which goes back to ideas of Risch, Norman, Davenport, Trager and others is a heuristic one and develops the integrand f over the differential field K in terms of multivariate rational functions such that

[pic]

Here v is a quotient of multivariate polynomials, the ui are logarithmic derivatives. Hence, the parallel integration method consists in making educated guesses for the ui and the denominator of v, as well as for the degree of its numerator; thus the calculation of algebraic quantities is reduced to the problem of solving a linear system of equations in order to determine the ci's and constant coefficients of the numerator of v. If the linear equations for the unknown coefficients has a solution, then an integral of f is found. Yet, the nonexistence of a solution does not always imply that f does not have an elementary integral over K but could mean that the guess was wrong. This makes the method although heuristic nevertheless very successful in some CAS such as Maple or Mathematica.

The original version of the program pmint (which stands for Poor Man's Integrator ) is a very short one with less than 100 lines of Maple code for integrating transcendental elementary or special functions.

Because powerful Maple functions (e.g. convert) used do not have a direct correspondence in Mathematica the task rewriting pmint into Mathematica code is not straightforward. However, the Mathematica version of pmint was formulated in 2007 by Luis A. Medina and Alexander Pavlyk at Wolfram Research Inc. Many examples which have been tested with the Maple version of pmint are investigated with the Mathematica implementation of pmint .

It turns out that the Poor Man's Integrator is able to compute even integrals which cannot be handled by the standard integrator implemented in Mathematica or Maple. This is due to the fact that pmint is applicable to special functions (such as Airy, Bessel, Lambert W, Wright Omega and Whittaker functions etc. ). The discussion of the current paper shows the strengths but also some drawbacks of this heuristic approach which can be considered as an extension of the well-established Risch algorithm

Method of Parallel Integration

In May 2005 Manuel Bronstein / INRIA presented on the 9th Workshop on Computer Algebra at LIT/JINR in Dubna his Maple version of the so-called "Poor Man's Integrator" a heuristic approach.

In his monography on "Symbolic Integration I" [13], published in 2005, there is a new chapt. 10 devoted to "Parallel Integration" - a method which goes back to ideas of R. Risch, A.C.Norman, J.H. Davenport, B.M. Trager and others; for references see [1-12]. This method attempts to handle all the generators of the differential field containing the integrand "in parallel", i.e. all at once rather than considering only the topmost one. Hence, this method has been called either the "new Risch algorithm" [5], the "Risch-Norman algorithm" [7,8,12] or the "parallel Risch algorithm" [9,10,11]; it is in fact heuristic rather than algorithmic as it may fail both in theory and practice in computing elementary antiderivatives.

As to [13] the general idea behind parallel integration is to avoid the recursive nature of the integration algorithm by viewing the differential field Κ containing the integrand f as a field of multivariate rational functions over its constants. The integrand f ∈ Κ belongs to a differential field of the form Κ = C(t1,... t n) where C = Const(Κ) and each ti is a transcendental function over C(t1,... t i-1).

As a consequence of the strong version of Liouville's Theorem (Bronstein 2005, Chap. 5.5.3, p.145 [13]), if the integral of f is elementary over Κ, then there are

v ∈ Κ, with c1,... cm algebraic numbers over C and u1,... um ∈ Κ(c1,... cm )* such that

[pic] (1)

Since v is a quotient of multivariate polynomials in t1,... tn and the ui's are assumed (without loss of generality) to be polynomials in t1,... tn (logarithmic derivative identity), the parallel method consists in making educated guesses for the ui 's and the denominator of v, as well as for the degree of its numerator. This reduces the problem of solving equation (1) to finding the c i 's and the constant coefficients of the numerator of v, which can be achieved with elementary linear algebra.

If the linear equations for the unknown constants c i have a solution, then an integral of f is found. However, the nonexistence of a solution does not always imply that f does not have an elementary integral over Κ; it could just mean that the guess was wrong. All parallel approaches published so far share the property that there are functions which have elementary antiderivatives that cannot be found with that variant. As a consequence, parallel integration is a convenient heuristic method which is significantly easier to implement than a complete integration algorithm. As will be shown below it is very successfully implemented for CAS such as Maple and Mathematica.

Mathematica code for pmint

The Maple implementation of the heuristic parallel integration, pmint, was published by M. Bronstein in May 2005 ("") and contains less than 100 lines of code. pmint which stands for "Poor Man's Integrator" is a program for integrating transcendental elementary or special functions

The analogous Mathematica implementation is more involved due to the fact that some powerfull built-in Maple procedures (such as convert etc. ) will have no direct correspondence in Mathematica. Yet, in 2007 at Wolfram Research Inc. (WRI) a Ph.D. student Luis A. Medina and Alexander Pavlyk [14] have formulated the Mathematica version of pmint, thus credit should be given to them. Various examples have been tested by the author under the most recent Mathematica V6.0.2

Subsequently, the Mathematica code for the procedures below is listed and - where available - compared with the corresponding Maple code :

(1) Code : pmint, pmIntegrate, tryIntegral, getSpecial, myFactors, enumerateMonoms, splitFactor, deflation

--------------------------------------------- pmint ------------------------------------------------------

Here is the Mathematica module :

Clear[pmint];

pmint[f_,x_,opt:OptionsPattern[{Extension→0}]]:=

Module[{ff,si,li,lin,lout,q,d,l,lv,ld,ls,fint,lc,t,terms,list},

ff=Together[TrigToTan[f]];

list=ToIndetsAndDerivation[ff,x,t];

If[list?$Failed,Return[INT[f]]];

{ff,terms,lv,ld}=list;

q=denD[ld];

ld=q*ld;

ls=DeleteCases[Map[getSpecial[#,Thread[terms→Rest[lv]]]&,terms],Null];

ToTerms[ pmIntegrate[ff,lv,ld,q,ls,OptionValue[Extension]],

terms,Rest[lv]]/.{tan→Tan}

]

For comparison the corresponding Maple procedures will be shown.

pmint := proc(f,x)

local ff, si, li, lin, lout, ld, q, d, l, vars, dx, ls, fint, lc;

ff := eval(convert(f, tan)); # convert trigs to tan

si := select(proc(d) diff(d,x) 0 end, indets(ff));

si := select(proc(d) diff(d,x) 0 end, indets(map(diff, si, x) )) union si;

li := [op(si)]; # list of terms in integrand and its derivative

lin := [seq(d=`tools/genglobal`(x), d=li)]; # substitution terms->

indets

lout := [seq(rhs(d)=lhs(d), d=lin)]; # substitution indets->terms

ld := subs(lin, map(diff, li, x)); # derivatives of all the terms

q := lcm(seq(denom(d), d=ld)); # denominator of the total derivation

l := [seq(normal(q * d), d=ld)]; # q * derivatives of all the terms

vars := map(lhs, lout);

dx := totalDerivation(vars, l); # vector field dx = q * d/dx

ls := [seq(getSpecial(d, lin), d=li)]; # list of known Darboux for dx

fint := subs(lout, pmIntegrate(subs(lin, ff), dx, q, vars, ls));

lc := select(proc(d) convert(d,string)[1]="_"; end, indets(fint,

name));

subs({seq(d = 0, d=lc minus si)}, fint);

end;

---------------------------------------------------------- pmIntegrate --------------------------------------------

Mathematica module :

Clear[pmIntegrate];

pmIntegrate[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 pmIntegrate."];

splq=splitFactor[q,lv,ld];

s=First[splq];

Scan[If[Last[#],s=s*First[#]]&,ls];

x=First[lv];

df=Denominator[f];

spl=splitFactor[df,lv,ld];

cden=s*First[spl]*deflation[Last[spl],lv,ld];

dg=1+totalDeg[s,lv]+

Max[ totalDeg[Numerator[f],lv],totalDeg[Denominator[f],lv]];

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=tryIntegral[ f,lv,ld,q,cand,lunk,

First[spl], (* normal part of df *)

Last[spl], (* special part of df *)

Last[splq], (* special part of q *)

ls,ext];

If[ First[sol],INT[f],Last[sol]]

]

Maple procedure

pmIntegrate := proc(f, d, q, vars)

local ls, splq, s, ff, df, spl, cden, dg, monomials, cand, lunk, sol, i;

if nargs = 5 then ls := args[5]; else ls := []; fi;

splq := splitFactor(q, d);

s := splq[1];

for i to nops(ls) do if ls[i][2] then s := s*ls[i][1]; fi; od;

ff := normal(f); df := denom(ff); spl := splitFactor(df, d);

cden := s * spl[1] * deflation(spl[2], d);

dg := 1 + degree(s) + max(degree(numer(ff)), degree(denom(ff)));

monomials := [op(enumerateMonoms(vars, dg))];

cand := add('_A'[i] * monomials[i], i = 1..nops(monomials)) / cden;

lunk := { seq('_A'[i], i = 1..nops(monomials)) };

sol:= tryIntegral(f, d, q, vars, cand, lunk,

spl[1], spl[2], splq[1], ls, 0);

if sol[1] then sol := tryIntegral(f, d, q, vars, cand, lunk,

spl[1], spl[2], splq[1], ls, I); fi;

if sol[1] then Int(f); else sol[2]; fi;

end;

---------------------------------------------------------- tryIntegral --------------------------------------------

Mathematica module :

Clear[tryIntegral];

tryIntegral[f_,lv_List,ld_List,q_,cand_,lunk_,l1_,l2_,l3_,ls_,k_]:=

(* In this application the "k" is a variable or variables to adjoint *)

Module[{candlog,p,candidate,i,sol,DD,BB,lu,dens},

DRPrint["Visit tryIntegral."];

DD=totalDerivation[lv,ld];

candlog=Union[myFactors[l1,k],myFactors[l2,k],

myFactors[l3,k],First /@ ls];

candidate=cand + Total[Table[BB[i] Log[candlog[[i]]],

{i,Length[candlog]}]];

lu=Union[lunk,Table[BB[i],{i,Length[candlog]}]];

DRPrint[" Applying 'funky way' of Together..."];

sol=totalDerivation[lv,Together[ld/q]];

{sol,dens}=Reap[Collect[f-sol[candidate], Alternatives @@ lu,

With[{tmp=Together[#]},Sow[Denominator[tmp]]; tmp]&]

];

dens=PolynomialLCM @@ Flatten[dens];

DRPrint[" Applying Collect..."];

sol=Collect[sol,Alternatives@@lu,Together[dens*#]&];

sol=DeleteCases[PolyCoeffs[sol,lv],0];

DRPrint[" Applying Outer..."];

dens=Outer[D,sol,lu];

DRPrint[" Solving linear equation with a ",

First[Dimensions[dens]]," x ", Last[Dimensions[dens]]," matrix."];

sol=Quiet[LinearSolve[dens,-sol/.Thread[lu→0]]];

If[Head[sol]===LinearSolve,sol={}];

If[sol=!={},sol=Thread[lu→sol]];

{sol==={},(candidate/.sol)/.Thread[lu→0]}

]

Maple procedure.

tryIntegral := proc(f, d, q, vars, cand, lunk, l1, l2, l3, ls, K)

local candlog, p, candidate, i, sol;

candlog := [op({myfactors(l1, K), myfactors(l2, K), myfactors(l3, K)}

union { seq(p[1], p=ls) })];

candidate := cand + add('_B'[i] * log(candlog[i]),

i = 1..nops(candlog));

sol := solve({coeffs(numer(normal(f - d(candidate)/q)), {op(vars)})},

lunk union { seq('_B'[i], i = 1..nops(candlog)) });

[evalb(sol=NULL), subs(sol,candidate)];

end;

---------------------------------------------------------- getSpecial --------------------------------------------

Mathematica module :

(* To get the known Darboux polynomials *)

getSpecial[f_tan,ru_]:={1+(f/.ru)^2,False};

getSpecial[f_Tanh,ru_]:= With[{ff=f/.ru},Sequence @@ {{1+ff,False},

{1-ff,False}}];

getSpecial[f_ProductLog,ru_]:={f/.ru,True};

getSpecial[___]:=Null;

Maple procedure.

getSpecial := proc(f, l) local p; # return known Darboux polynomials

p := op(0,f);

if p = `tan` then [1+subs(l,f)^2, false];

elif p = `tanh` then [1 + subs(l,f), false], [1 - subs(l,f), false];

elif p = `LambertW` then [subs(l,f), true];

else NULL; fi;

end;

---------------------------------------------------------- myfactors --------------------------------------------

Mathematica module :

myFactors[p_,Algebraic]:=

Module[{f,var},

( var=Variables[#];

If[Length[var]?1,

Apply[Sequence,( First[var]-(First[var]/.

{ToRules[ Roots[#?0,First[var]] ]})

)],#] )& /@ myFactors[p,I]

];

myFactors[p_/;!NumericQ[p],0]:= Drop[First/@ FactorList[p],1];

myFactors[p_?NumericQ,_]:={p};

myFactors[p_/;!NumericQ[p],extension_]:=

Drop[First/@ FactorList[p,Extension→ extension],1];

Maple procedure.

myfactors := proc(p, K) local l, fact;

if K = 0 then l := factors(p); else l := factors(p, K); fi;

seq(fact[1], fact=l[2]);

end;

--------------------------------------------- enumerateMonoms --------------------------------------------------

Mathematica module :

(* Calculate all monomials in variables "lv" such that TOTAL deg is ≤ "deg" *)

enumerateMonoms[{},deg_]:={1};

enumerateMonoms[lv_List,deg_]:=

Module[{i,v=Most[lv]},

Union[enumerateMonoms[v,deg],

Sequence @@ Table[Last[lv]^i*enumerateMonoms[v,deg-i],{i,deg}]]

]

Maple procedure.

enumerateMonoms := proc(vars, d) local n, x, i, v, s, w;

n := nops(vars);

if n = 0 then {1}; else

x := vars[n];

v := [seq(vars[i], i = 1..n-1)];

s := enumerateMonoms(v, d);

for i to d do s := s union {seq(x^i*w,w=

enumerateMonoms(v,d-i))};

od; s;

fi;

--------------------------------------------- splitFactor --------------------------------------------------------

Mathematica module :

splitFactor[p_,lv_List,ld_List]:=

Module[{theT,c,hn,hs,S,q,qn,qs,DD},

theT=mainVar[p,lv];

If[theT===$Failed,Return[{p,1}]];

DD=totalDerivation[lv,ld];

{c,q}=PolyContentPP[p,theT];

{hn,hs}=splitFactor[c,lv,ld];

S=PolynomialQuotient[PolynomialGCD[q,DD[q]],

PolynomialGCD[q,D[q,theT]],theT];

If[Exponent[S,theT]?0,Return[{hn p,hs}]];

{qn,qs}=splitFactor[PolynomialQuotient[q,S,theT],lv,ld];

{hn qn,S hs qs}

]

Maple procedure.

splitFactor := proc(p, d) local si, x, c, q, spl, s, splh;

si := select(proc(z) d(z) 0 end, indets(p,name));

if si = {} then RETURN([1,p]) fi;

x := si[1];

c := content(p, x, 'q');

spl := splitFactor(c, d);

s := normal(gcd(q, d(q)) / gcd(q, diff(q, x)));

if degree(s) = 0 then RETURN([spl[1], q * spl[2]]); fi;

splh := splitFactor(normal(q / s), d);

[spl[1] * splh[1] * s, spl[2] * splh[2]];

end;

--------------------------------------------- deflation -----------------------------------------------------------

Mathematica module :

(* Calculate the deflation of p with respect to the "new derivation". *)

deflation[p_,lv_List,ld_List]:=

Module[{theT,c,q,DD},

theT=mainVar[p,lv];

If[theT===$Failed,Return[p]];

DD=totalDerivation[lv,ld];

{c,q}=PolyContentPP[p,theT];

deflation[c,lv,ld]*PolynomialGCD[q,DD[q,theT]]

]

Maple procedure.

deflation := proc(p, d) local si, x, c, q;

si := select(proc(z) d(z) 0 end, indets(p,name));

if si = {} then RETURN(p) fi;

x := si[1];

c := content(p, x, 'q');

deflation(c, d) * gcd(q, diff(q, x));

end;

(2) Conversion of trig functions to tan : TrigToTan, Derivative

--------------------------------------------- TrigToTan ------------------------------------------------------

Mathematica module :

Clear[TrigToTan,tan,exp];

TrigToTan[ee_]:=

Module[{ToTan,e=ee/.Sec[x_]^2→1+Tan[x]^2},

ToTan[Sin[x_]]=With[{y=Together[x/2]},2 tan[y]/(1+tan[y]^2)];

ToTan[Cos[x_]]=With[{y=Together[x/2]},(1-tan[y]^2)/(1+tan[y]^2)];

ToTan[Sec[x_]]=With[{y=Together[x/2]},(1+tan[y]^2)/(1-tan[y]^2)];

ToTan[Sec[x_]^2]:=1+tan[x]^2;

ToTan[Csc[x_]]=With[{y=Together[x/2]},(1+tan[y]^2)/(2 tan[y])];

(* ToTan[Tan[x_]]=With[{y = Together[x/2]},2 tan[y]/(1-tan[y]^2)]; *)

ToTan[Tan[x_]]:=tan[x];

ToTan[Cot[x_]]=With[{y=Together[x/2]},(1-tan[y]^2)/(2 tan[y])];

e//.{x:(_Sin|_Cos|_Tan|_Cot|_Sec|_Csc)ƒToTan[x]}

]

In Maple there is a powerful built-in procedure convert(expr,form) which transforms expr into another form.

The function convert is used to convert an expression from one form to another. Some of the conversions are data-type conversions, others are form or function conversions. For the latter one a set of optional arguments are given to perform the conversion in different manners, e.g. convert(f,tan)which converts trigonometric functions into expressions involving only the function tan.

--------------------------------------------- Derivative ----------------------------------------------------------

The Mathematica module :

Derivative[1][tan]=Function[1+tan[#]^2];

Derivative[1][exp]=exp;

exp[a_+b_]:=exp[a] exp[b];

exp[Log[x_]]:=x;

is an extension of the first derivatives of tan and exp with additional rules for the exponential function

(3) Other auxillary functions are : totalDerivation, PolyContentPP, PolyCoeffs, mainVar, myVariables, RationalFunctionQ, td, totalDeg, SubSet.

--------------------------------------------- totalDerivation -------------------------------------------------

Mathematica module :

totalDerivation[lv_List,ld_List]:=

Function[u,ld.Table[D[u,k],{k,lv}],Listable]

Maple procedure.

totalDerivation := proc(lv, ld)

proc(f) local fp, i;

fp := 0; for i to nops(ld) do fp := fp + ld[i] * diff(f, lv[i]); od;

fp;

end;

end;

--------------------------------------------- PolyContentPP ----------------------------------------------------

Mathematica module :

(* Returns the content and the primitive part of A as a polynomial in t *)

PolyContentPP[A_,t_]:=

With[{q=FactorTermsList[A,t]},

{ Times @@ Most[q],Last[q] }

]

--------------------------------------------- PolyCoeffs ----------------------------------------------------------

PolyCoeffs[p_,vars_]:=

Block[{t},

Union @ Flatten[ CoefficientArrays[p,vars] /.

t_SparseArray :> ArrayRules[t] [[All,2]] ]

]

--------------------------------------------- mainVar ----------------------------------------------------------

mainVar[p_,lv_List]:=

Catch[ Scan[ If[Exponent[p,#1]>0,Throw[#1]]&,Reverse[lv]];

$Failed

]

--------------------------------------------- myVariables --------------------------------------------------------

myVariables[e:(_Plus|_Times)]:=

Union[Flatten[myVariables/@ Apply[List,e]]];

myVariables[HoldPattern[Power][e_,a_Integer]]:= myVariables[e];

myVariables[e_/;NumericQ[e]]:= {};

myVariables[HoldPattern[Power][x_,p_Plus]]:=

Union[Flatten[myVariables[x^Apply[List,p]]]];

myVariables[e_]:= {e};

SetAttributes[myVariables,{Listable}];

--------------------------------------------- RationalFunctionsQ ------------------------------------------

RationalFunctionQ[expr:(_Plus|_Times),x_]:=

Fold[Function[{pr,el},pr && RationalFunctionQ[el,x]],True,

List @@ expr];

RationalFunctionQ[HoldPattern[Power][e_,a_Integer],x_]:=

RationalFunctionQ[e,x];

RationalFunctionQ[x_,x_]:= True;

RationalFunctionQ[expr:Except[_Times|_Plus],x_List]:=

Fold[#1 && (expr === #2||FreeQ[expr,#2])&,True,x];

RationalFunctionQ[expr:Except[_Times|_Plus],x:Except[_List]]:=

FreeQ[expr,x];

--------------------------------------------- td ------------------------------------------------------------------------

td[mono_]:= Total[Exponent[mono,myVariables[mono]]];

td[mono_,vars_]:= Total[Exponent[mono,vars]];

--------------------------------------------- totalDeg -------------------------------------------------------------

totalDeg[poly_Plus]:= Max[td /@ List @@ poly];

totalDeg[poly_]:= td[poly];

totalDeg[poly_Plus,vars_]:= Max[td[#,vars]& /@ List @@ poly];

totalDeg[poly_,vars_]:= td[poly,vars];

--------------------------------------------- SubSet -----------------------------------------------------------------

SubSet[A_,B_]:=

Module[{a=Tally[A][[All,1]],b=Tally[B][[All,1]]},

Length[Intersection[a,b]]? Length[a]

];

(4) Building monomials extension : ToIndetsAndDerivation, ToTerms, denD

--------------------------------------------- ToIndetsAndDerivation --------------------------------------

Clear[ToIndetsAndDerivation];

ToIndetsAndDerivation[f_,x_,t_]:=

Module[{vars,terms,deriv,newF,limitList,list,F,dterms,count=0},

F=f//.{Exp→exp};

terms=Select[myVariables[F],PolynomialQ[#,x]?False &];

dterms=Union[Select[ Flatten[myVariables[D[terms,x]]],

(PolynomialQ[#,x] ? False &)],terms];

While[ Not[SubSet[dterms,terms]] && count≤15,

terms = Union[terms,dterms];

dterms=Union[Select[Flatten[myVariables[D[terms,x]]],

(PolynomialQ[#,x] ? False &)],terms];

count=count+1

];

If[count>15, Return[$Failed] ];

deriv=D[terms,x];

vars=Map[ t,Range[Length[terms]] ];

limitList=Thread[terms -> vars];

newF=F /. limitList;

deriv= deriv /. limitList /.

{HoldPattern[Power][a_,b_Plus] :>

Times @@ ((a^Apply[List,b]) /. limitList)};

If[ True || Fold[#1 && RationalFunctionQ[#2,Join[{x},vars] ]&,

True,deriv],

{newF,terms,Join[{x},vars],Join[{1},deriv]}, $Failed]

(* Else *)

]//.{exp[y_] :> Exp[y]};

Here f is the integrand, x the variable of integration and t a dummy variable whose purpose is to name the new variables as t1, t2 ,… . The output of the procedure ToIndetsAndDerivation is :

• the integrand f in terms of new variables,

• to adjoin the variables in order to calculate the integral. This information is needed when changing back the result in terms of the original variables.

• The last two elements of the output are of the form {x, t1, t2],…} and {1,D[t1], D[t2],…}. These elements are needed to obtain totalDerivation.

--------------------------------------------- ToTerms ---------------------------------------------------------------

ToTerms[F_,terms_List,vars_List]:= F/.Thread[vars -> terms]

It is assumed that "derivations" are the list of derivatives of the ti 's. This list should be given at the time, when denD is called. In particular, it takes the last element of ToIndetsAndDerivation[f,x,t].

--------------------------------------------- denD --------------------------------------------------------------------

denD[derivations_List]:= PolynomialLCM @@ Denominator[Together[derivations]];

Finally, some auxillary definitions are needed for simplification of the expressions etc.

(5) Transformations : LambertRule,log2ArctanRule,pureFctConversion

Because Lambert's W function is implemented in Mathematica as the function ProductLog[z] a substitution rule is introduced.

--------------------------------------------- LambertRule --------------------------------------------------------

LambertRule = {ProductLog[z_] -> lambertW[z]};

The inverse of the function f(w) = w eW is called the Omega function. The Wright Omega function is defined in terms of the Lambert W function with argument ©x .

--------------------------------------------- Omega function --------------------------------------------------

Ω[x_]:=LambertW[©x]

In order to recast the logarithmic representation of the arctan function into the usual representation a sequence of replacement rules has to be defined

--------------------------------------------- log2ArctanRule --------------------------------------------------

log2ArctanRule [result_]:=

((( ( result/.{Log[a_]-Log[b_] -> Log[a /b]} )/.

{Log[(-™+z_)/ ™+z_] -> Log[- ((1+™ z)/(1-™ z))]} )/.

{Log[-c_] -> Log[c]+Log[-1]} )/.

{Log[(1+™ z_)/(1-™ z_)] -> 2™ ArcTan[z]}

)//Simplify

The following rule achieves the conversion between the two forms of Mathematica's pure function representation with Function[par,body] and body[#]& .

--------------------------------------------- pureFctConversion ---------------------------------------------

pureFctConversion =

Literal[ Function[ par_, body_ ] ] /; FreeQ[ body, #] :>

Unevaluated[ Function[ body ] /. par → # ];

and finally a debugging utility

--------------------------------------------- DRPrint ---------------------------------------------------------------

Clear[DRPrint];

DRPrint[w___]/;($DebugRisch===True):= Print[w];

Demonstration of pmint

In order to demonstrate the potential of the heuristic parallel integrator pmint a selection of examples will given. There are several categories such as:

a) Elementary functions

b) Rational functions

c) Trigonometric functions

d) Log-Exp functions

e) Liouvillian special functions ( Error function)

f) Airy functions

g) Bessel functions

h) Lambert W function

i) Wright Omega function

As will be shown below some examples involving Airy, Bessel, Lambert W, Wright Omega and Whittaker W functions where pmint gives a result cannot be dealed with either by Maple's built-in integrator int or/and by Mathematica's built-in procedure Integrate. Where available the results of the built-in integrators obtained by Maple and/or Mathematica will be compared with the result of pmint.

(a) Elementary functions

In the case of elementary functions such as (non-rational) polynomials involving powers of x or combinations of trigonometric functions pmint reproduces the same results as found with the built-in integrators of Mathematica or Maple.

Generally, there arise some problems to cast the results in the same form. Thus, as shown by the subsequent example it needs some additional efforts, i.e. a sequence of Mathematica commands, until the result looks the same.

res1 = pmint[x2 Csc[x3] Sec[x3],x]

=> -(1/3) Log[-1+Tan[x3/2]]+1/3 Log[Tan[x3/2]]-1/3 Log[1+Tan[x3/2]]

The result which looks quite complicated can be simplified by a sequence of transformation rules

res11=(((res1//TrigFactor )//.{Log[a_]-Log[b_]→Log[a /b]}

//FullSimplify)/.{Log[(-(1/2))*a_]→-Log[-2]+Log[a]}//Expand);

res11//.{(c=Select[res11,FreeQ[#,x]&])→ C1}

=> 1/3 Log[Tan[x3]]+C1

and agrees with the result of Integrate

res2 = Integrate[x2 Csc[x3] Sec[x3],x] //Simplify

up to an complex integration constant C1=-(1/3) (Log[2]+ ™ π ) as can be easily proofed by differentiating both results D[res1-res2,x] //Simplify => 0 .

(b) Rational functions

Handling rational polynomial functions seems to be no problem for pmint

[pic]

=> [pic]

However, the obviously simple integration [pic] with result ArcTan[x] only works with extension to complex numbers ™ by the additional option of pmint , i.e. Extension→Algebraic.

[pic]

=> -(1/2) ™ (Log[-™+x]-Log[™+x])

With a sequence of replacement rules log2ArctanRule the logarithmic terms are finally casted into the function ArcTan.

res1 //log2ArctanRule

=> ArcTan[x]+ π/2

which is the result of Integrate up to an integration constant [pic] .

(c) Trigonometric functions

Moreover, pmint can handle integrands which are rational polynomials of trigonometric functions and powers of x.

However, pmint cannot cope with the following simple example, because [pic] has no representation in terms of transcendentals t i. The result of Integrate is an elliptic integral of 2nd kind [pic]. But, if the series expansion of the integrand, i.e. series[[pic],{x,0,14}]//Normal, is used pmint will deal with the integrand up to arbitrary order. The result in terms of a power series

[pic]

coincides with the series expansion of [pic].

Straightforward calculcation of the integral [pic] fails because it requires an extention to rationals for factorization. Therefore, with the additional option Extension ( Algebraic which is passed to FactorList the procedure pmint will find a solution

res1 = pmint[1/(2+Sin[x]),x,Extension→Algebraic]

=> [pic]

Alternatively, the Extension ( (-1)1/3 could have been chosen but this choice would require some kind of preknowledge of the result; therefore the option Extension ( Algebraic is more general. In order to obtain the result in the same form given by the built-in integrator first the roots {(-1)1/3,(-1)2/3} have to be casted into [pic] by means of the substitution rule algNumbRule :

algNumbRule = Thread[{(-1)1/3,(-1)2/3} ( ComplexExpand[{(-1)1/3,(-1)2/3}] ];

res1a =(res1/.algNumbRule//FullSimplify)/.{Log[arg_](Log[2 arg]}//Simplify

With the help of the subsequent rewrite rule

(((res1a /. {1+2Tan[x/2] ( [pic]z,-1-2Tan[x/2] ( -[pic]z}//Simplify )/.

{Log[a_]-Log[b_] ( Log[a/b]})/.

{Log[(™-z)/( ™+z)] ( -(Log[1-™ z]-Log[1+™ z])}//ExpToTrig )/.

{z ( (1+2Tan[x/2])/[pic]}

the result finally agrees with that of Integrate.

=> [pic]

Another intricate example is

[pic]

A straightforward calculation fails again but the calculation is doable if the integrand is casted into a rational polynomial using the substitution Tan[x/2] → z. Because of this substitution the transformed integrand fct(z) has to be divided by an additional factor [pic] . Hence using Apart and ComplexExpand the integrand becomes

[pic] and [pic]

Now,

ϑ1 = pmint[f1,z,Extension → Algebraic]//. {z → Tan[x/2]} //Simplify

=> 1+Cos[x]-2™ Log[-™+Tan[x/2]]+2™ Log[™+Tan[x/2]]+Log[1+Tan[x/2]]+Sin[x]

The result of the second part ϑ2 is quite lengthy

ϑ2 = ((pmint[f2,z,Extension → Algebraic] //Simplify )//. {z → Tan[x/2]})

ϑ2//Short[#,4]&

=> (22/3 (?1?) Log[-1+1/6 (1+™[pic]) (54-6 [pic])1/3+((1-™[pic]) (9+?1?)1/3)/62/3+Tan[x/2]])/(-264 21/3 35/6 [pic]+?4?+(9-[pic])2/3 (-35/6 (13 ™+3 [pic]) (2 (9+[pic]))1/3+(39+3 ™[pic]) ?1?+6 (9+3 ™[pic]+™ [pic]+[pic])))+?1?/(31/6 (?1?))+((?1?) ?1?)/(6?1? (?1?))+(4 (?1?))/?1?

but it is an explicit solution. Yet, interesting is the comparison with Integrate. The straightforward application of Integrate leads to a result much shorter which, however, is given as an implicit representation in terms of Mathematica's RootSum objects

Κ = Integrate[(Sin[x] Cos[x] -Sin[x]2)/(1+Sin[x]+Cot[x]) ,x]//FullSimplify

=> 2 x+Cos[x]+Log[Sec[x/4]2]+Log[Cos[x/2]+Sin[x/2]] +Sin[x] +

8 RootSum[-7+296 #1+704 #12+5632 #13&,

Log[-25+2 (9+88 #1 (-1+8 #1)) Tan[x/4]+25 Tan[x/4]2] #1& ]

Numerical approximations of the results both of Κ and Τ1+Τ2 agree up to integration constants.

(d) Log-Exp functions

Here is an example where pmint gives a solution

[pic]

=> Log[©2 x+[pic] Log[-1+x]3]

but Integrate fails.

(e) Liouvillian special functions

Liouvillian special functions are rational polynomials containing powers of the error functions erf(x) together with [pic] :

[pic]

=> [pic]

The result of integration is the same for pmint and the integrators built into Maple and Mathematica. However, for a slightly more complicated case where the numerator [pic] is replaced by [pic] neither pmint nor the implemented integrators of Maple or Mathematica can deal with this case.

(f) Airy functions

For the subsequent examples which involve AiryAi and AiryAiPrime functions pmint is able to find solutions while Integrate fails in all cases.

pmint[x AiryAiPrime[x] Cos[AiryAi[x]]+Sin[AiryAi[x]],x]

=> [pic]

Here is another integral involving combinations of AiryAi, AiryAiPrime, Log and ArcTan :

[pic]//

Simplify

=> AiryAi[x] ArcTan[x] Log[x]

Another integral which pmint can solve but not Integrate is

[pic]

=> [pic]

Finally

[pic]

=> [pic]

In order to obtain this result the CPU time (with respect to the computer being used) is quite long.

(g) Bessel functions

A straightforward integration for the quotient of Bessel functions [pic] fails both for pmint and Integrate.

[pic]

The reason why pmint cannot directly cope with this integrand and hangs up is that every differentiation introduces new orders of BesselJ functions in this specific case

NestList[ D[#,x]&,BesselJ[n,x],4] //Simplify//

TraditionalForm //ColumnForm

=> Jn(x)

1/2 (Jn-1(x)-Jn+1(x))

1/4 (Jn-2(x)-2 Jn(x)+Jn+2(x))

1/8 (Jn-3(x)-3 Jn-1(x)+3 Jn+1(x)-Jn+3(x))

1/16 (Jn-4(x)-4 Jn-2(x)+6 Jn(x)-4 Jn+2(x)+Jn+4(x))

In order to circumvent this feature of generating additional orders of BesselJ functions one applies instead a recurrence relation for Bessel functions according to [14] to obtain a simple differentiation rule.

( D[{BesselJ[n,x],BesselJ[n+1,x]},x]//

Simplify`RecurrenceShift[#,TargetFunctions→ BesselJ[n,x]]& //

Simplify /. {BesselJ[n,x]→t1, BesselJ[n+1,x]→t2}

=> [pic]

This result serves as a guidance for the definition of t1' and t2' in terms of pure functions .

{Derivative[1][t1] = Function[z,[pic]]/. pureFctConversion ,

Derivative[1][t2] = Function[z,[pic]]/. pureFctConversion}

=> { [pic] , [pic]}

The great advantage of pmint is that it suffices to know the derivative of a function (which needs to be in the monomial extension of transcendental functions) only in order to be able to perform the integration. It is assumed that the differential is closed on that extension. Therefore, because the derivative of

D[n Log[x]-Log[t1[x]],x]//Simplify

=> [pic]

turns out to be [pic] the procedure pmint is able to calculate the integral.

[pic]

=> [pic]

The same strategy can also be used for the calculation of

[pic]

where a direct calculation again fails but with the help of the auxillary functions t1 and t2 and their 1st derivatives given above)the integration turns out to be

[pic]

=> [pic]

(h) Lambert W function

In Mathematica the Lambert W function is identical with the function ProductLog. Hence, either LambertW or ProductLog can be used.

The function w = ProductLog[z] is the principal solution for w of z = w ©w. (The kth solution is denoted by ProductLog[k,z]). Moreover, ProductLog[z] satisfies the differential equation [pic].

Both pmint and Integrate can handle the following integral :

pmint[Sin[ProductLog[x]],x]//FullSimplify

=> [pic]

For the next example

[pic]

a straightforward integration with pmint does not work. Because pmint has the default OptionsPattern[{Extension ( 0}] this example fails due to the fact that certain polynomials cannot be factorized. However, it suffices to enlarge the extension to Extension ( i in order to obtain the solution.

[pic]

=> (™/2)(Log[-™+x ProductLog[x]]-Log[™+x ProductLog[x]])

With the help of the transformation rule log2ArctanRule the logarithms are finally simplified to

=> (/2+ArcTan[x ProductLog[x]]

For the following integral a solution only exists with pmint whereas Integrate does not.

pmint[(2 LambertW[x2] Cos[LambertW[x2]] (a x+LambertW[x2])+

a x (1+LambertW[x2] )+2 LambertW[x2])/((1+LambertW[x2])

(a x+LambertW[x2]) x),x] // FullSimplify

=> Log[a x+ProductLog[x2]]+Sin[ProductLog[x2]]

(i) Wright Omega-function

The Wright Omega function is defined in terms of the LambertW function with argument ©x. Therefore, we define

Ω[x_]:=LambertW[©x]

Again, pmint gives a solution whereas Integrate fails :

pmint[Ω[x],x]//Simplify

=> -1+ProductLog[©x]+1/2 ProductLog[©x]2

This applies to the following example too :

[pic]

=> Log[x+ProductLog[©x]]+Sin[ProductLog[©x]]

Again, the built-in integrator Integrate cannot not find a solution whereas pmint does.

Conclusion :

It has to be emphasized that the Poor Man's Integrator ( pmint ) is not intended to be as complete as the large integrators based on the Risch algorithm which are built into Maple and Mathematica. The very small size of pmint makes it easy to port this procedure to any computer algebra system or library capable of factoring multivariate polynomials. Because it is applicable to special functions (such as Airy, Bessel, Lambert W and Wright Omega functions etc. ), the Poor Man's Integrator is capable to compute integrals some of which cannot be handled by existing integrators. Therefore, pmint is not intended to be a replacement for existing integrators - rather an extension, or a cheap and powerful alternative for any computer algebra project.

References :

[1] R. Risch (1968 ) : On the integration of elementary functions which are built up using algebraic operations.

Research Report SP-2801/002/00, System Development Corp., Santa Monica, CA, USA

[2] R. Risch (1969) : Further results of elementary functions.

Research Report RC-2402, IBM Research, Yorktown Heights, NY, USA

[3] R. Risch (1969) : The problem of integration in finite terms.

Transactions of the American Mathematical Soc., 139 : 167-189

[4] R. Risch (1970) : The solution of problem of integration in finite terms.

Bulletin of the American Mathematical Soc., 76 : 605-608

[5] A.C.Norman, P.M.A.Moore (1977) : Implementing the new Risch Integration Algorithm,

Proceedings of the 4th International Colloquium on Advanced Computing Methods in Theoretical Physics, 99-110.

[6] Risch (1979) : Algebraic properties of the elementary functions of analysis.

American Journal of Mathematics, 101: 743-759

[7] S.H.Harrington (1979) : A new symbolic integration system in reduce,

The Computer Journal 22, 127-131.

[8] J.Fitch (1981) : User-based integration software,

Proceedings of SYMSAC'81, ACM Press, 245-248.

[9] J.H.Davenport (1982) : On the Parallel Risch Algorithm (I),

Proceedings of EUROCAM'82, LNCS 144, Springer, 144-157.

[10] J.H.Davenport (1982) : On the Parallel Risch Algorithm (III): Use of Tangents,

SIGSAM Bulletin 16, 3-6.

[11] J.H.Davenport, B.M.Trager (1985) : On the Parallel Risch Algorithm (II),

ACM Transactions on Mathematical Software 11, 356-362.

[12} K.Geddes, L.Stefanus (1989) : On the Risch-Norman Integration Method and its Implementation in Maple, Proceedings of ISSAC'89, ACM Press, 212-217.

[13] M. Bronstein (2005) : Symbolic Integration I, Transcendental Functions (2nd Ed.) , Springer-Verlag Berlin Heidelberg

[14] A. Pavlyk / WRI (Nov. 2007) : Private communications

................
................

In order to avoid copyright disputes, this page is only a partial summary.

Google Online Preview   Download