maxima-admin@math.utexas.edu wrote:
> I implemented the method of M. Bronstein & S. Lafaille,
> for solving second-order linear DEs; for a description
> of this method, see their ISSAC'2002 article
>
> http://www-sop.inria.fr/cafe/Manuel.Bronstein/bronstein-eng.html
>
> 0. It might not be the fastest method, but the method of
> M. Bronstein & S. Lafaille is the most general method I
> know of for solving second-order linear DEs in terms of
> hypergeometric functions. To speed up solving easy cases,
> I added a few solvers that work by pattern matching. If
> you know of something better, let me know.
>
> 1. If you have done similar work, you might let me know.
> My code allows new methods to be added without fuss. If
> we could agree on a standard for calling these methods,
> our work could be blended.
>
> 2. My code checks each solution. When the solution
> is say y = Bessel_j(0,x), the checking code needs to
> apply the Bessel function recursion relation to crunch
> the result to zero. Sigh. I "fixed" this by changing
> the gradef for the Bessel functions. Now,
> diff(Bessel_j(n,x),x) --> dBessel_j(n,x) and
> diff(dBessel_j(n,x) --> <junk> * Bessel(n,x).
> This isn’t good; I really need something else…
> (But I think each solution should be checked.)
>
Barton,
I faced a similar issue when checking Bessel function solutions
of Riccati equations. The code is in share/contrib/diffequations.
1. I wrote a simple ode checking function ode_check() in
file contrib_ode.
2. There is a test harness in the tests subdirectory.
3. There are tests that return Bessel function solutions in
file rtest_ode1_riccati.mac. I defined identities
/* Need these indentities to check solutions */
besjident(n,x):=bessel_j(n,x)=2*(n-1)*bessel_j(n-1,x)/x-bessel_j(n-2,x);
besjident(n,x):=bessel_j(n,x)=2*(n-1)*bessel_j(n-1,x)/x-bessel_j(n-2,x);
besiident(n,x):=bessel_i(n,x)=-2*(n-1)*bessel_i(n-1,x)/x+bessel_i(n-2,x);
besiident(n,x):=bessel_i(n,x)=-2*(n-1)*bessel_i(n-1,x)/x+bessel_i(n-2,x);
beskident(n,x):=bessel_k(n,x)=2*(n-1)*bessel_k(n-1,x)/x+bessel_k(n-2,x);
beskident(n,x):=bessel_k(n,x)=2*(n-1)*bessel_k(n-1,x)/x+bessel_k(n-2,x);
then use them as required
/* Bessel function solution - non-integer order b*c<0 */
ans:ode1_riccati_original(b:1,c:-1,m:2,y,x);
y=((bessel_j(-5/4,x^2/2)*%c+bessel_j(-3/4,x^2/2))*x^2+bessel_j(-1/4,x^2/2)*%c)/((bessel_j(-1/4,x^2/2)*%c+bessel_j(1/4,x^2/2))*x);
ode_check(eqn,ans);
(((bessel_j(-1/4,x^2/2)+bessel_j(-9/4,x^2/2))*%c+bessel_j(1/4,x^2/2)+bessel_j(-7/4,x^2/2))*x^2+5*bessel_j(-5/4,x^2/2)*%c+3*bessel_j(-3/4,x^2/2))/(bessel_j(-1/4,x^2/2)*%c+bessel_j(1/4,x^2/2));
expand(ev(%,besjident(1/4,x^2/2),besjident(-1/4,x^2/2)));
0;