New package for drawing direction fields with gnuplot

Dnia Fri, 19 Nov 2010 15:11:14 -0500, Mark H Weaver napisa?(a):

> Adam Majewski <adammaj1 at o2.pl> writes:
>> soln_at draw field lines, so how can I draw equipotential lines ?
> 
> With the attached patch to drawdf.mac, does the following command do
> what you want?
> 
> drawdf(x*x+y*y, [x,-2,2], [y,-2,2],
>   field_color=gray, proportional_axes='xy, aspect_ratio=1, color=red,
>   solns_at(makelist([k/2,-k/2],k,-3,3)), color=green,
>   equipots_at(makelist([k/2,k/2],k,-3,3)))$
> 
> Note that equipots_at uses the same numerical integrator as solns_at, so
> you may sometimes have to mess around with options such as duration and
> tstep to get good results.
> 
>     Mark
> 
> 
> --- drawdf-orig.mac	2010-11-17 20:09:09.000000000 -0500 +++ drawdf.mac
> 2010-11-19 14:57:11.000000000 -0500 @@ -204,6 +204,13 @@
>    if curve=[] then curve:[pt],
>    [ 'point_type='none, funmake('soln_curve, [curve, 'soln_arrows =
>    soln_arrows]) ] )$
>  
> +/*
> + * Ideally, this should probably detect and remove overlaps in the
> curve. + */
> +equipotential_plot(pt) := block(
> +  [ dxdt:dydt, dydt:-dxdt, soln_arrows:false, soln_direction:'both ], +
>  solution_plot(pt) )$
> +
>  linearize_at(pt, derivs, vars) := block(
>    [ subs: map("=", vars, pt) ],
>    subst(subs, apply('matrix, outermap('diff, derivs, vars))))$
> @@ -348,6 +355,10 @@
>      then pre_graphics: append(reverse(solution_plot(args(param))),
>      pre_graphics) else if operatorp(param, 'solns_at)
>      then for pt in list_of_points_from(param) do pre_graphics:
>      append(reverse(solution_plot(pt)), pre_graphics)
> +    else if operatorp(param, 'equipot_at) +    then pre_graphics:
> append(reverse(equipotential_plot(args(param))), pre_graphics) +    else
> if operatorp(param, 'equipots_at) +    then for pt in
> list_of_points_from(param) do pre_graphics:
> append(reverse(equipotential_plot(pt)), pre_graphics)
>      else if operatorp(param, 'saddle_at)
>      then pre_graphics: append(reverse(saddle_plot(args(param))),
>      pre_graphics) else if operatorp(param, 'saddles_at)
> _______________________________________________ Maxima mailing list
> Maxima at math.utexas.edu
> http://www.math.utexas.edu/mailman/listinfo/maxima

I have done :

 patch < drawdf-equipot.diff
(Stripping trailing CRs from patch.)
patching file drawdf.mac
Hunk #1 succeeded at 206 (offset 2 lines).
Hunk #2 succeeded at 357 (offset 2 lines).


Yes, it works fine.


I have also tried :

f(zx,zy):=block([i:0,zx2,zy2,cx:0,cy:0,ER2:4,iMax:20],
zx2:zx*zx,zy2:zy*zy,
while zx2+zy2<ER2 and i<iMax do 
(zy:2*zx*zy+cy,zx:zx2-zy2+cx,zx2:zx*zx,zy2:zy*zy,i:i+1),
return(i));


drawdf(f(x,y), [x,-2,2], [y,-2,2]);

Unable to evaluate predicate y^2+x^2<4
#0: f(zx=x,zy=y)
 -- an error. To debug this try: debugmode(true);

Adam