(*

## Mathematica Package for Imitating Calculator-Generated Slope Fields

*)

`(* :Title: Slope Field Plots for First Order ODEs in Two Variables *)(* :Author: Jack Courtney*)(*:Summary:        This package plots slope fields in the Cartesian plane for first order        ODEs of the form y'[x]==f[x,y], where f is a function of x and y.*)(* :Context: Graphics` *)(* :Package Version: 1.0 *)(* :Mathematica Version: 2.2.3 *)(* :Copyright: Copyright (c) 1995, 1996 by Jack Courtney.  All Rights Reserved. *)(* :Name: Graphics`SlopeField` *)(* :Keywords:        slope field, differential equation, ode*)(* :Requirements: None. *)(* :Warnings: None. *)(* :Sources: *)BeginPackage["Graphics`SlopeField`"]PlotSlopeField::usage="PlotSlopeField[ derivative_, xvals_L\ist, yvals_List, opts___ ], where derivative is the name of \the function of x and y whose value is y'[x], and the list\s xvals and yvals are fully qualified ranges. For example, \if y'[x_]:=h[x,y[x]], then the field of slopes in |x|<1, |y\|<2 is produced by SlopeField[ h, {-1,1,0.1}, {-2,2,0.2} ].  \The (optional) third argument in each of the two ranges g\ives the sample spacing, and opts is an (optional) list of \Graphics options.";PlotSlopeFieldAndSamples::usage="PlotSlopeFieldAndSamples[ \derivative_, xvals_List, yvals_List, throughPts_List, opts_\__ ] calls SlopeField[ derivative_, xvals_List, yvals_List, \opts___ ]; then, for each of the points listed in the vari\able list 'throughPts', plots the exact solution of y'[x]==\derivative[x,y[x]] that contains that point.  For example, \if throughPts={{1,2},{1,3}}, then the graphs of the exact s\olutions that contain (1,2) and (1,3) are both added to the \graph of the slope field.";Begin["`Private`"]Off[ General::spell1 ];ct=Context[];ToCleanString[ exp_ ]:=   StringReplace[ ToString[ exp ], ct->"" ];PlotSlopeField[ deriv_, xvals_List, yvals_List, opts___ ]:=  Module[ { LineThicknessRatio=0.001             (* i.e., 1/1000 of the width of the plot *)          , DefaultSeparationRatio=0.04          (* i.e., 1/25 of the width of the plot   *)          , ArrowLengthRatio=0.2                 (* i.e., 1/5 of the cell dimension       *)          , dx,dy, xinc,yinc          , lines, label, df, r          }        , If[ Length[xvals]<3            , xinc=N[ (xvals[[2]]-xvals[[1]])*DefaultSeparationRatio ]            , xinc=N[ xvals[[3]] ]            ];          If[ Length[yvals]<3            , yinc=N[ (yvals[[2]]-yvals[[1]])*DefaultSeparationRatio*(AspectRatio/.Options[Graphics]) ]            , yinc=N[ yvals[[3]] ]            ];          dx=xinc*ArrowLengthRatio //N; (* the width of a grid cell *)          dy=yinc*ArrowLengthRatio //N; (* the height of a grid cell *)          r=dy/dx;                      (* the aspect ratio of a grid cell *)          (*          ** Each arrow is centered in a grid cell, with both head and tail lying on the boundary.          *)          x=xvals[[1]]+xinc/2 //N;          lines={};          While[ x<N[xvals[[2]]]               , y=yvals[[1]]+yinc/2 //N;                 While[ y<N[yvals[[2]]]                      , If[ Abs[df=deriv[x,y]]>r                          , AppendTo[ lines, {{x-dy/df,y-dy},{x+dy/df,y+dy}} ]                          , AppendTo[ lines, {{x-dx,y-dx*df},{x+dx,y+dx*df}} ]                          ];                        y+=yinc                      ];                 x+=xinc               ];          Off[ Unset::norep ];          y=.; x=.; y'[x]=.; y[x]=.;          On[ Unset::norep ];          label=ToString[StringJoin[ "Slope Field for "                                   , ToCleanString[ y'[x_]==deriv[x,y[x]] ]                                   , ",\n \n "                                   ]];          Show[ Graphics[ {Thickness[LineThicknessRatio], Map[ Line,lines,{1} ]}                        , Axes->True                        , AxesLabel->{"x","y"}                        , PlotLabel->label                        , PlotRange->{Take[xvals,2],Take[yvals,2]}						, DefaultFont->{"Courier",8}                        ]]        ];PlotSlopeFieldAndSamples[ deriv_, xvals_List, yvals_List, throughPts_List, opts___ ]:=  Module[ { IC, soln, label,plot0,plot1}        , soln=DSolve[ y'[x]==deriv[x,y[x]], y[x], x][[1,1,2]];          IC=Apply[(soln/.C[1]->Solve[(soln/.x->#1)==#2,C[1]][[1,1,2]]//N)&                  , throughPts                  , {1}                  ]/.x->s;          label=ToString[StringJoin[ "Slope Field for "                                   , ToCleanString[ y'[x_]==deriv[x,y[x]] ]                                   , ",\n       With Special Solutions\n "                                   ]];          Block[ {\$DisplayFunction=Identity},                 plot0=PlotSlopeField[ deriv, xvals, yvals ];                         plot1=Plot[ Evaluate[IC/.s->x]  (* Show the slope field with          *)                           , {x,xvals[[1]],xvals[[2]]}      (*   the sample solns. plotted on it. *)                           , PlotRange->{Take[xvals,2],Take[yvals,2]}                           ]               ];          Show[ plot0, plot1              , DefaultFont->{"Courier",8}              , PlotLabel->label              ]        ];Protect[ PlotSlopeField, PlotSlopeFieldAndSamples ];On[ General::spell1 ];End[];EndPackage[];(*(*Demo1:*)Off[ General::spell ];Off[ General::spell1 ];yPrime[x_,y_]:=-3*x*y;                                   (* Define the diffeq y'[x_] = g[x,y]          *)PlotExactSolnsThruPts={{0,2},{0,-1}};                    (* Plot exact solutions that pass through ... *)PlotXRange={-3,3,0.3};PlotYRange={-2,2,0.3};plot=PlotSlopeField[yPrime,PlotXRange,PlotYRange];SetOptions[\$Output,PageWidth->80];Print[ "USAGE:  ", OutputForm[PlotSlopeField::usage] ];  (* Print out the caption for the plot; *)SetOptions[\$Output,PageWidth->51];                       (*   i.e., the usage of SlopeField     *)On[ General::spell1 ];On[ General::spell ];*)(*(*Demo2:*)Off[ General::spell ];Off[ General::spell1 ];yPrime[x_,y_]:=-3*x*y;                                   (* Define the diffeq y'[x_] = g[x,y]          *)PlotExactSolnsThruPts={{0,2},{0,-1}};                    (* Plot exact solutions that pass through ... *)PlotXRange={-3,3,0.3};PlotYRange={-2,2,0.3};PlotSlopeFieldAndSamples[yPrime,PlotXRange,PlotYRange,PlotExactSolnsThruPts];SetOptions[\$Output,PageWidth->80];Print[ "USAGE:  ", OutputForm[PlotSlopeFieldAndSamples::usage] ];SetOptions[\$Output,PageWidth->51];On[ General::spell1 ];On[ General::spell ];*)`
`(9/12/96) jac`