Plotting direction fields for ODEs



Hi,
I have finally finished the implementation of a plotdf function for
Maxima. You can see the details in http://fisica.fe.up.pt/maxima/

I now that this is not a good time to incorporate new functionalities,
but in this case it implies just a few lines of lisp, to recover some
functions that were already present in the tcl/tk scripts but that
were not being used.

If you want to, I can send some patches now, for the current CVS
version, or I can wait until further notice. If Jim grants me access
(villate@sourceforge) I could also commit the patches by myself.

See the documentation for the new function at the end (PLOT_OPTIONS
keeps growing and I could now incorporate a lot more options for plot2d
and openplot_curves :)

Cheers,
Jaime
----
(%i1) ? plotdf;


Info from file /usr/local/info/maxima.info:
 - Function: PLOTDF (expr,...,options,..)
 - Function: PLOTDF ([expr1,expr2],...,options,..)
     Plots the direction field for the first order differential equation
                 dy
                 -- = EXPR
                 dx
     where EXPR is an expression depending on x, y and a set of
     parameters that must be given numerical values with the PARAMETERS
     option. The second form plots the direction field for the
     autonomous system:
                 dx             dy
                 -- = EXPR1     -- = EXPR2
                 dt             dt
     where EXPR1 and EXPR2 can only depend on variables x and y, apart
     from a list of given numerical parameters.

     Several options can be given in the command, or in the menu that
     will appear when you click on the upper left corner of the plot
     window.  Integral curves can be obtained by cliking on the plot,
     or with the option TRAJECTORY_AT. The direction of the integration
     can be control with the DIRECTION option. The number of
     integration steps is given by NSTEPS and the time interval between
     them is TSTEP. The Adams Moulton method is used for the
     integration; it is also posible to switch to an adaptive
     Runge-Kutta 4th order method.

     The menu in the plot window has the following options: Zoom, will
     change the behavior of the mouse so that it allows you to zoom in
     on a region of the plot. Each click near a point magnifies the
     plot, keeping the center at the point you clicked. Depressing the
     SHIFT key while clicking zooms in the opposite direction. To
     resume computing trajectories click on Integrate. To change the
     differential equation, or change various settings, click on Config
     and enter new values in the entry windows, and then click on
     Replot in the main menu bar.

     Holding the right mouse button down allows you to drag (translate)
     the plot sideways or up and down. Additional parameters such as
     the number of steps (nsteps), the initial t value (tinitial), and
     the x and y centers and radii, may be set under the  Config menu.

     You may print to a postscript printer, or save the plot as a
     postscript file, by clicking on save. To change between printing
     and saving see the Print Options under Config.

     To obtain the direction field for the equation diff(y,x) = x - y²,
     and the solution with initial condition y(-1) = 3, we can use
     the command:

          plotdf(x-y^2,[xfun,"sqrt(x);-sqrt(x)"],[trajectory_at,-1,3],[direction,forward],[yradius,5],[xcenter,6]);
     plots the direction field for the equation diff(y,x) = x - y^2,
     and the solution with initial condition y(-1) = 3.  The graph also
     shows the function y = sqrt(x).

          plotdf([y,-k*x/m],[parameters,"m=2,k=2"],[sliders,"m=1:5"],[trajectory_at,6,0]);
     Will show the direction field of a harmonic oscillator, and
     the integral curve through (x,y) = (6,0), with a slider that
     will allow you to change the value of m interactively (k is fixed
     at 2).

          plotdf([y,-(k*x + c*y + l*x^3)/m],[parameters,"k=-1,m=1.0,c=0,l=1"],[sliders,"k=-2:2,m=-1:1"],[bbox,-3,-3,3,3],[tstep,0.1]);
     will show the direction field of the Duffing equation.

     See PLOT_OPTIONS for more examples.