(*
*)
(* :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