Initial- and boundary-value problems in Maxima.

>Jaime's code in dynamics.mac looks rather nice. If the changes
>needed to make it "decent" are to add adaptive step size and
>error control, it would be great to add that.
>
>If the term "decent" in Gregor Berman's original post
>also means "much faster" then a different tactic may be necessary.
>Typically this means writing out the function-evaluation parts into
>fortran-syntax files, compiling them and linking with fortran
>programs, and communicating with Maxima / plotting/ etc.  This
>kind of setup has been used in the past, but tends to be more
>delicate -- for example, not all of the operating systems/lisps
>that support Maxima will manage such connections to fortran in
>the same way.  Also, the programs that spew out fortran syntax
>may not be as sophisticated as needed.
>
>An intermediate possibility exists, which is to use lisp
>rather than Maxima language / [even compiled Maxima language]
>which should be considerably faster.
>
>There are other points to note:  the arithmetic in Maxima can be
>altered to be higher precision, and so the treatment of sensitive
>parts of the solver can be modified to take smaller steps and also
>do the arithmetic to (arbitrarily) high precision if this would
>help.
>
>RJF

English is not my native language, I am really sorry if the term "decent"
was not polite or appropriate. It didn't meant to be rude or anything. Now,
rk() is nice, but implements a rather introductory algorithm, still useful
for some experiments (students often ask for such a neat Runge-Kutta method
implementation, because it is what they learn in their first lesson in
Numerical Analysis).
However, constant step size is never used for serious work, for two reasons.
First, it is not possible to know in advance how many steps will be needed
for an accurate solution; that always depends from the specific problem at
hand. Second, in most cases, one needs to use a very small constant step to
get an accurate solution, although that small step is only needed for a
small part of the integration range, while a much larger step is more than
sufficient for the rest. This makes constant step methods very slow.

I understand that the main focus of Maxima has been symbolic computations
rather than numerical methods. However, it is not easy to find a realistic
problem, taken from Physics or other sciences, for which a complete analytic
solution is possible. Even a very simple mechanical system, such as the
pendulum, is described mathematically in a way that numerical solutions are
necessary. For that reason, nearly all Computer Algebra Systems have
built-in support for numerical computations, using implementations of
numerical algorithms in the background. Most of them are translations of
well-tested Fortran programs, mainly posted in http://www.netlib.org.

I understand how rk() works. I believe a small correction should be made at
line 13 of the code (line 391 of dynamics.mac), which reads

if (not(numberp(last(apply(f_rk_4,cons(t0,xxx)))))) then

This line catches an error in the initial conditions only if it happens at
the last initial condition in list. For example, initial conditions like
[1,0,k] will be detected as erroneous (provided k is not defined), but
initial conditions like [1,k,0] will not.  I think this line should be
replaced with

if member(false,map(numberp,apply(f_rk_4,cons(t0,xxx)))) then

This ensures any non-numerical initial condition will be detected.

I'll be back with a Maxima implementation of a initial-value problems
solver, with built-in step size control soon.