Is there a numerical ODE solver facility in Maxima

I wanted to say thank you to everyone, particularly Leo Butler, who
embarrassed me by demonstrating that this is a problem I am familiar with,
and Jaime Villate, who showed me how to use rk and with this problem, and
gave me a fairly in depth example of how to use it.

I ended up using rk in a shooting method in the problem I posted and a more
difficult one.  Both of these problems are from liquid crystal physics.
This method worked as well as could be expected, I suppose.  This is my
first time using this method, and it turns out have to be kind of smart
about how you do it.  Several times I ran into machine precission problems
when dealing with machine sized floats (of course I could switch to bfloats,
but it turned out that I didn't care enough to make this work in every
corner case).  By the way, is there any documentation on the bfloat
implementation?

Would it be helpful to have a shooting method ODE solver in Maxima, I think
maybe so for people like me.  I might be able to generalize my code into
something that can be used for other problems, but I will finally need to
learn how to write "good" Maxima code.

In my class, students were strongly encouraged to use Mathematica and
examples and solutions were provided in Mathematica.  It was hard to keep up
with the other students who use Mathematica, as you might expect.  Some of
this was due to the (assumedly) extensive work that Wolfram's programers
have put into ironing out corner cases in various procedures.  I.e. where
the other students just type NDSolve[blah], I actually had to think about
how to solve the problem.  I am glad that I now know how to solve this kind
of problem now, even without Maxima.  A typical physicist, I would argue,
could care less about the route to the answer as long as he can be
reasonable sure it is correct, or at least that he is using the "industry
standard" for solving differential equations.  I am not saying this out of
some sense of need to tell people what way to steer Maxima, I am just saying
that this is my experience.

Zach


On Sun, Apr 5, 2009 at 3:08 PM, Jaime Villate <villate at fe.up.pt> wrote:

> On Sex, 2009-04-03 at 15:24 -0600, Zach wrote:
> > Here is the equation:
> >
> > diff( f(x), x, 2) = -k * sin( f(x) ) * cos( f(x) ) /* k is just some
> > positive constant */
> >
>
> > rk(['diff(g, x), -1*sin(f)*cos(f)], [f,g], [0,0]/*initial value*/, [x,
> > -.5, .5, .01]); /*for k=1*/
>
> No, that's not how you should use rk.
>
> First, notice that the state variables are f, and its derivative, g:
>  (%i1) vars: [f, g]$
>
> The derivatives of those two state variables (for k=1) are:
>  (%i2) derivs: [g, -sin(f)*cos(f)]$
>
> They define the components of the phase velocity. Plot the phase
> portrait using:
>  (%i3) plotdf(derivs,vars,[f,-5,5],[g,-3,3]);
>
> The result (see figure attached) resembles a pendulum.
> There are equilibrium points at (f=n*%pi/2, g=0), n=0, +-1, +-2,...
> We can see in the phase portrait that the equilibrium points at m*%pi
> are centers (stable) and those at (2*m+1)*%pi/2 are saddle points
> (unstable).
>
> You can also deduce the nature of the equilibrium points by looking at
> the eigenvalues of the jacobian matrix:
>  (%i3) J: jacobian (derivs, vars);
>  (%i4) eigenvalues (ev (J, f=0, g=0));
>  (%o4) [[- %i, %i], [1, 1]]
> (purely imaginary eigenvalues proves that (0,0) is a center).
>  (%i5) eigenvalues(ev(J,f=%pi/2,g=0));
>  (%o5) [[- 1, 1], [1, 1]]
> (two real eigenvalues with the same absolute value mean that (%pi/2,0)
> is a saddle point).
>
> The curves that join the saddle points in the phase portrait constitute
> what is called a heteroclinic cycle. If you give initial conditions
> within those regions, the solution will be an oscillation around the
> stable equilibrium point.
>
> It is not worth giving initial conditions [0,0] for [f.g], since [0,0]
> is an equilibrium point. If you want to study an oscillatory solution,
> use, for instance:
>  (%i6) sol1: rk(derivs, vars, [0, 0.5], [x, 0, 20, 0.1])$
> To see the function f(x):
>  (%i7) plot2d([discrete, makelist([sol1[i][1], sol1[i][2]],
>            i, 1, length(sol1))], [xlabel,"x"], [ylabel, "f(x)"])$
>
> And a solution outside of the region of oscillatory solutions:
>  (%i22) sol2: rk(derivs, vars, [0, 1.1],[x, 0, 20, 0.1])$
>  (%i23) plot2d([discrete, makelist( [sol2[i][1], sol2[i][2]],
>        i, 1, length(sol1))], [xlabel,"x"], [ylabel, "f(x)"])$
>
>
> Regards,
> Jaime
>
>