Is there a numerical ODE solver facility in Maxima

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

-------------- next part --------------
A non-text attachment was scrubbed...
Name: portrait.pdf
Type: application/pdf
Size: 24353 bytes
Desc: not available
Url : http://www.math.utexas.edu/pipermail/maxima/attachments/20090405/d5269024/attachment-0001.pdf