I've written a couple of functions for plotting vector fields in Maxima. My
Maxima skill is only modest, so I suspect experts may find more elegant ways
to code this. The labeling of the axes, if nothing else, could be improved,
since my code labels them with lisp names (i.e. $X instead of x).
If you try this code and find the arrow heads on your 3D plots look wrong, I
believe it is due to a bug in Maxima, which I reported (with a patch) today:
https://sourceforge.net/tracker/?func=detail&atid=104933&aid=3136608&group_id=4933
You should be able to get acceptable plots if you start Maxima fresh,
and
try your 3D plots right away, before doing any other plotting. (In truth, I
suspect there might be a second bug here--not respecting any head_length
commands until after draw2d has run--but it works.)
/* plot_vector_field( F, X, Y, ... )
* Draws a 2D vector field
* F - a vector containing field components
* X - name and bounds of first coordinate
* Y - name and bounds of second coordinate
* optional parameters:
* scale=1 auto-scaling of arrows (default)
* scale=0 no auto-scaling
* scale=* adjust arrows shorter -- numbers between 0 and 1
* or longer -- numbers greater than 1
* other parameters pass through to draw environment
* Example:
* load(draw)$
* plot_vector_field( [-y,x], [x,-1,1],[y,-1,1] )$
* plot_vector_field( [-y,x], [x,-1,1],[y,-1,1], scale=0 )$
*/
plot_vector_field( F, X, Y, [args] ) := block(
[vars, G,P,Q, points, grid_d,u,scale, range,xrange,yrange, v],
scale : assoc('scale,args,1), args:delete( 'scale=scale, args ),
/* create function versions of each componont and of the field itself */
vars : [ X[1],Y[1] ],
P : G( vars, ev(F[1])), P : subst( lambda, G, P ),
Q : G( vars, ev(F[2])), Q : subst( lambda, G, Q ),
G : lambda( [z], [apply(P,z), apply(Q,z)] ),
/* create a list of points to base arrows at */
points : listify( cartesian_product(
setify( makelist( X[2] + (X[3]-X[2])*i/10, i, 0, 10 )),
setify( makelist( Y[2] + (Y[3]-Y[2])*i/10, i, 0, 10 )) )),
/* diagonal length of a grid square */
grid_d : sqrt((X[3]-X[2])^2 + (Y[3]-Y[2])^2)/10.0,
/* u is the divisor that shortens arrows to fit in grid squares */
u : max( apply( max, makelist( abs(P(k[1],k[2])), k, points )) /
(X[3]-X[2]) * 10,
apply( max, makelist( abs(Q(k[1],k[2])), k, points )) /
(Y[3]-Y[2]) * 10 ),
if scale#0 and u>grid_d/50.0 then u:scale/u else u:1,
/* xrange and yrange need to be large enough to contain head
* and tail of each arrow (or some won't appear in output) */
range : lambda( [z], [ apply(min,z), apply(max,z) ] ),
xrange : range( flatten(makelist( [k[1],k[1]+P(k[1],k[2])], k, points ))),
yrange : range( flatten(makelist( [k[2],k[2]+Q(k[1],k[2])], k, points ))),
/* generate the vectors to draw and set an appropriate arrow
* head length for each one */
v : flatten( makelist( [head_length=u*sqrt(G(k).G(k))/5+grid_d/100,
vector(k, u*G(k) )], k, points )),
/* draw! */
apply( draw2d,
append([color=blue],[line_width=1],[head_type='nofilled],[head_angle=20],['xrange=xrange],['yrange=yrange],[xlabel=X[1]],[ylabel=Y[1]],
args, v ))
)$
/* plot_vector_field3d( F, X, Y, Z, ... )
* Draws a 3D vector field
* F - a vector containing field components
* X - name and bounds of first coordinate
* Y - name and bounds of second coordinate
* Z - name and bounds of third coordinate
* optional parameters:
* scale=1 auto-scaling of arrows (default)
* scale=0 no auto-scaling
* scale=* adjust arrows shorter -- numbers between 0 and 1
* or longer -- numbers greater than 1
* other parameters pass through to draw environment
* Example:
* load(draw)$
* plot_vector_field3d( [-y,x,z], [x,-1,1],[y,-1,1],[z,-1,1] )$
* plot_vector_field3d( [-y,x,z], [x,-1,1],[y,-1,1],[z,-1,1], scale=0 )$
*/
plot_vector_field3d( F, X, Y, Z, [args] ) := block(
[vars, G,P,Q,R, points, grid_d,u,scale, range,xrange,yrange,zrange, v],
scale : assoc('scale,args,1), args:delete( 'scale=scale, args ),
/* create function versions of each componont and of the field itself */
vars : [ X[1],Y[1],Z[1] ],
P : G( vars, ev(F[1])), P : subst( lambda, G, P ),
Q : G( vars, ev(F[2])), Q : subst( lambda, G, Q ),
R : G( vars, ev(F[3])), R : subst( lambda, G, R ),
G : lambda( [z], [apply(P,z), apply(Q,z), apply(R,z)] ),
/* create a list of points to base arrows at */
points : listify( cartesian_product(
setify( makelist( X[2] + (X[3]-X[2])*i/5, i, 0, 5 )),
setify( makelist( Y[2] + (Y[3]-Y[2])*i/5, i, 0, 5 )),
setify( makelist( Z[2] + (Z[3]-Z[2])*i/5, i, 0, 5 )) )),
/* the diagonal length of a grid square */
grid_d : sqrt((X[3]-X[2])^2 + (Y[3]-Y[2])^2 + (Z[3]-Z[2])^2)/5.0,
/* u is the divisor that shortens arrows to fit in grid squares */
u : max( apply( max, makelist( abs(P(k[1],k[2],k[3])), k, points )) /
(X[3]-X[2]) * 5,
apply( max, makelist( abs(Q(k[1],k[2],k[3])), k, points )) /
(Y[3]-Y[2]) * 5,
apply( max, makelist( abs(R(k[1],k[2],k[3])), k, points )) /
(Z[3]-Z[2]) * 5 ),
if scale#0 and u>grid_d/50.0 then u:scale/u else u:1,
/* xrange,yrange,zrange need to be large enough to contain head
* and tail of each arrow (or some won't appear in output) */
range : lambda( [z], [ apply(min,z), apply(max,z) ] ),
xrange : range( flatten(makelist( [k[1],k[1]+P(k[1],k[2],k[3])], k, points
))),
yrange : range( flatten(makelist( [k[2],k[2]+Q(k[1],k[2],k[3])], k, points
))),
zrange : range( flatten(makelist( [k[3],k[3]+R(k[1],k[2],k[3])], k, points
))),
/* generate the vectors to draw and set an appropriate arrow
* head length for each one */
v : flatten(makelist( [head_length=u*sqrt(G(k).G(k))/10+grid_d/100,
vector(k, u*G(k))], k, points )),
/* draw! */
apply( draw3d,
append([color=blue],[line_width=1],[head_type='nofilled],[head_angle=10],['xrange=xrange],['yrange=yrange],['zrange=zrange],[xlabel=X[1]],[ylabel=Y[1]],[zlabel=Z[1]],
args, v ))
)$
I would like to contribute this code to the Maxima project to be licensed
under the same terms as Maxima itself.
Donald J Bindner