> From: Robert Dodier
>
>
> On 2/26/07, Billinghurst, David (RTATECH)
> <David.Billinghurst at riotinto.com> wrote:
>
> > (%i4)
> besselcrunch((-(1+2*n)*bessel_j(n,x))/x+bessel_j(1+n,x)+bessel
> _j(n-1,x));
> >
> > Maxima encountered a Lisp error:
> >
> > Expected a maxima function designator but got NIL.
>
> OK, I get the same error with the current version of Maxima
> (5.11.0cvs)
> and share/contrib/opsubst.lisp (r1.5).
>
> Going back through the revisions of opsubst, the same error happens
> with r1.4 but not with r1.3 (timestamp = 2006/05/06).
> However, with r1.3 the results I get are different from the ones shown
> in Barton's original message.
> (http://www.math.utexas.edu/pipermail/maxima/2006/001287.html)
>
> e.g.
> besselcrunch((-(1+2*n)*bessel_j(n,x))/x+bessel_j(1+n,x)+bessel
> _j(n-1,x));
> => -2*n*?%bessel_j(n,x)/x-?%bessel_j(n,x)/x+?%bessel_j(n+1,x)
> +?%bessel_j(n-1,x)
>
> (Uh-oh -- looks like there's a bug in MSIZE-ATOM or something
> ... rats.)
>
> It looks like the noun/verb distinction is causing the opsubst stuff
> to fail to recognize %BESSEL_J as the same as $BESSEL_J ...
> The following definition of OP-EQUALP (otherwise defined in
> src/compar.lisp)
> yields the expected result for the above example (namely 0).
>
> (defun op-equalp (e &rest op) (and (consp e) (consp (car e)) (some
> #'(lambda (s) (equal ($verbify (caar e)) s)) (mapcar #'$verbify
> op))))
>
> Dunno if that causes unexpected results elsewhere.
Thanks for reply. I too struck some problems with %BESSEL_J vs $BESSEL_J.
I have checked in some code to share/contrib/diffequations/contrib_ode.mac.
It is an amalgam of Barton's ideas and some code I had been tinkering with.
It works for all but two of the Bessel function solutions in the contrib_ode
testsuite.
/* This Bessel function simplification code can be improved,
but is does work fairly well on the contrib_ode testsuite.
Need to extend to cover cases like Kamke 1.14 and 1.24
*/
bessimp(e) := besksimp(besisimp(besysimp(besjsimp(e))))$
freeof_bessel(e):=freeof(nounify(bessel_j),nounify(bessel_y),
nounify(bessel_i),nounify(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 contains Bessel functions of order n-1 and n-2.
We separate the set of arguments into groups that have a common
artgument 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(nounify(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)
)
),
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)));
besksimp(ex):=besFsimp(ex,bessel_k,
lambda([n,x], 2*(n-1)*bessel_k(n-1,x)/x+bessel_k(n-2,x)));
NOTICE
This e-mail and any attachments are private and confidential and may contain privileged information. If you are not an authorised recipient, the copying or distribution of this e-mail and any attachments is prohibited and you must not read, print or act in reliance on this e-mail or attachments.
This notice should not be removed.