> # BESSEL FUNCTIONS

>

> # What does Maple know about Bessel functions?

> ?Bessel

> ?BesselZeros

>

> # The exercises occupying pp. 379-386 of D. Betounes,

> # Partial Differential Equations for Computational Science

> # (Springer, 1998)

> # show lots of things that computer algebra systems can do with Bessel functions.

>

> # For example, Ex. 5 challenges your CAS to evaluate the integrals that arise

> # as normalization constants in Fourier-Bessel series:

>

> Int(BesselJ(0, w*r)^2 * r, r);

> value(%);

>

> Int(BesselJ(3, w*r)^2 *r, r);

> value(%);

> Int(BesselJ(n, r)^2 *r, r);

> value(%);

> # Let's try to expand something in Bessel functions.

> h := x -> 1/(1 + x^2);

> # Use the notation for the 0th-order Fourier-Bessel series on p. 96 of notes.

> # Take r_0 = 1.

> c := k-> num(k)/den(k);

> den := k -> int( BesselJ(0, BesselJZeros(0,k)*r)^2 * r, r=0..1);

> den(1);

> evalf(%);

> # (It works as expected.)

> num := k -> int( BesselJ(0, BesselJZeros(0,k)*r) * h(r) * r, r=0..1);

> num(1);

> evalf(%);

> c(2);

> evalf(%);

> partialsum := (n,r) -> evalf(sum( c(k)*BesselJ(0,BesselJZeros(0,k)*r), k=1..n));

> plot(h, 0..1);

> plot(partialsum(1,r), r=0..1);

> plot({partialsum(2,r), h(r)}, r=0..1);

>

> plot({partialsum(6,r), h(r)}, r=0..1);

> plot({partialsum(20,r), h(r)}, r=0..1);

> # This looks much like a typical Fourier series, complete with Gibbs phenomenon

> # at the end where all the eigenfunctions vanish.

>

> # Let's try a different value of the Bessel index. This will cause trouble at

> # the origin, too.

>

> den := k -> int( BesselJ(3, BesselJZeros(3,k)*r)^2 * r, r=0..1);

> num := k -> int( BesselJ(3, BesselJZeros(3,k)*r) * h(r) * r, r=0..1);

> partialsum := (n,r) -> evalf(sum( c(k)*BesselJ(3,BesselJZeros(3,k)*r), k=1..n));

> plot({partialsum(1,r), h(r)}, r=0..1);

> plot({partialsum(10,r), h(r)}, r=0..1);

> plot({partialsum(50,r), h(r)}, r=0..1);

>