> From: Barton Willis
>
> A great project...Try something like this
>
> (1) use gatherargs to collect the arguments of the Bessel
> functions in an expression,
>
> (2) find the least (or greatest) Bessel function order,
>
> (3) use Bessel function recursion to collapse to reduce the Bessel
> function order
> (no more than two distinct orders that differ by an integer.
>
> It will be fun; I've made things like this before.
It was a great project ;-) Lots of fun with sets, partitions,
equivalence classes, etc; and then uglified by nasty little hack
to deal with bessel_k.
See bessimp (below) from share/contrib/diffequations/contrib_ode.mac).
This works for all the ode cases, but may not be robust enough for
everyday use.
There is also a function for simplifying expressions containing the
exponential integral expintegral_e using the recurrence (A&S 5.1.14).
expintegral_e(n+1,z)
= (1/n) * (exp(-z)-z*expintegral_e(n,z)) n = 1,2,3 ....
#######################################################################
/* This Bessel function simplification code can be improved,
but is does work well on the contrib_ode testsuite.
*/
bessimp(e) :=
if freeof_bessel(e) then
e
else
besksimp(besisimp(besysimp(besjsimp(e))))$
freeof_bessel(e):=freeof(bessel_j,bessel_y,bessel_i,bessel_k,e);
/* Simplify expressions containing bessel_F, where F in {j,y,i,k}.
Repeatedly find the highest order term F and substitute Fsub for F
When F is a Bessel function of order n and argument x, then
Fsub is contains Bessel functions of order n-1 and n-2.
We separate the set of arguments into classes that have identical
second arguments and orders that differ by integers, then work on
each
class until all the orders differ by less than 2.
*/
besFsimp(ex,F,Fsub) :=
block([arglist,ords,ord2,max,min,s,t,u,m,n,x,e:ex],
if freeof(F,ex) then return(ex),
argset:setify(gatherargs(ex,F)),
/* partition arguments into sets with same second argument and
orders that differ by an integer */
argset:equiv_classes(argset,lambda([x,y],
(second(x)=second(y)) and integerp(ratsimp(first(x)-first(y))))),
/* Iterate over each equivalance class */
for s in argset do (
x:second(first(s)), /* The common argument for s */
ords:map('first,s), /* List of orders */
t:first(ords), /* One element of the list of orders */
ord2:map(lambda([m],ratsimp(m-t)),ords), /* List of differences */
max:lmax(ord2), /* maximum order is t+max */
min:lmin(ord2), /* minimum order is t+min */
for m: max step -1 unless ratsimp(m-min)<2 do (
n:subset(ords,lambda([u],is(ratsimp(u-t-m)=0))),
n:if emptyp(n) then ratsimp(t+m) else first(n),
e:subst(apply(Fsub,[n,x]),apply(F,[n,x]),e),
e:ratsimp(e)
)
),
ratsimp(e)
)$
besjsimp(ex):=besFsimp(ex,bessel_j,
lambda([n,x], 2*(n-1)*bessel_j(n-1,x)/x-bessel_j(n-2,x)));
besysimp(ex):=besFsimp(ex,bessel_y,
lambda([n,x], 2*(n-1)*bessel_y(n-1,x)/x-bessel_y(n-2,x)));
besisimp(ex):=besFsimp(ex,bessel_i,
lambda([n,x],-2*(n-1)*bessel_i(n-1,x)/x+bessel_i(n-2,x)));
/* bessel_k requires special treatment due to the simplification
bessel_k(-n,x) => bessel_k(n,x), which upsets the recurrence
relationship. */
besksimp(ex):=block(
[bk:?gensym(),e],
e:opsubst(bk,bessel_k,ex),
e:besksimp_1(e,bk),
/* Recurrence relation for bessel_k */
e:besFsimp(e,bk,lambda([n,x], 2*(n-1)*bk(n-1,x)/x+bk(n-2,x))),
opsubst(bessel_k,bk,e)
);
/* This is a helper function for besksimp.
- the function bk has been substituted for bessel_k in ex
- collect all the arguments of each occurence of bk in e
- partition arguements into sets with same argument and where
either the orders m and n differ by integer or m and (-n)
differ by integer (so m+n is an integer).
- for the elements where m+n is an integer, substitute
bk(-n,x) for bk(n,x) in e
- return e
*/
besksimp_1(e,bk):=block(
[s,argset],
/* Break into equivalence classes */
argset:equiv_classes(setify(gatherargs(e,bk)),
lambda([x,y],
(second(x)=second(y))
and
(integerp(ratsimp(first(x)-first(y)))
or integerp(ratsimp(first(x)+first(y)))))),
/* For each equivalnce class */
for s in argset do block(
[x,ords,t,n],
x:second(first(s)), /* The common argument for s */
ords:map('first,s), /* Set of orders */
/* t is an element of ords, If ords are numbers then t is the
largest */
t:first(ords), if every(numberp,ords) then t:lmax(ords),
/* Find the subset that do not differ by an integer from t
and change sign of the order for those orders */
for n in subset(ords,lambda([i],integerp(ratsimp(t+i)))) do (
e:subst(bk(-n,x),bk(n,x),e) )
),
e
);