# # The user routines in this package are # # fieldplot: create a direction field plot for a first order ODE # phaseplot: create a phase plot with flow lines for two first order ODE's # firsteuler: solves a first order ODE using Eulers method # impeuler: solves a first order ODE using the "improved Euler method" # rungekutta: numerical solution to a system of first order ODEs using # a 4'th order Runge-Kutta method. The rungekuttahf routine # is a hardware floating point version of rungekutta # # There is a Maple V help file included for each of these routines. # After reading in the procedures, to get help for phaseplot type ?phaseplot # # June 1992. # Fixed arrow routines which would not scale arrows in extreme cases. # Changed (improved?) the shape of the arrows. # Added type checking to improve the error messages the user would # see when the parameters were incorrectly input, thanks to Michael # Monagan for his help on this. # # November 1991. # Added help files for rungekutta, rungekuttahf, firsteuler, impeuler. # Fixed a bug in rungekuttahf which was restricting accuracy to current # Digits setting, rather than the accuracy of evalhf. # # April 1991. # This file contains Maple procedures which work with Maple V. They are # used to analyze systems of first order differential equations. # # Hard copy documentation for the procedures is available as chapter 5 # of a book published by Brooks/Cole Publishing Company called Maple V # Flight Manual. # # Thanks to George Labahn of the Waterloo group for sharing his code, # parts of which helped to speed the following code. Also, the technique # of putting arrows on vectors is his. Michael Monagan has been very # helpful testing the routines and making suggestions for general # improvements. Others in the Waterloo group have also been helpful for # testing and suggestions. # # Daniel Schwalbe # 6635 County Road 101 # Corcoran, MN 55340 # # e-mail dschwalb@daisy.uwaterloo.ca # # `help/text/fieldplot` := TEXT( `FUNCTION: fieldplot - create a direction field with flow lines `, `for first order differential equation.`, ``, `CALLING SEQUENCE:`, ` fieldplot(f, h, v)`, ` fieldplot(f, h, v,{inits},...)`, ``, `PARAMETERS:`, ` f - function of two variables where D(y)(t)=f(t,y).`, ` h - horizontal range`, ` v - vertical range`, ` {inits} - set of initial values for flow lines.`, ``, `SYNOPSIS: - A typical call to the fieldplot function is`, `fieldplot(f,a..b,c..d), where f is a real function of two`, `variables and a..b,c..d specifies the horizontal and vertical real`, `range. For each point, (tp,yp), of a grid a short line segment is`, `drawn with slope f(tp,yp).`, ``, `- Each initial point is specified as a list, [t0,y0]. A fourth`, `order Runge/Kutta scheme is used to calculate the solution to`, `D(y)(t)=f(t,y) through the point (t0,y0) over the entire`, `horizontal range. In particular, the flow line will be`, `erratic over any interval where a solution does not exist.In`, `this case there is an option to use an adaptive numerical scheme`, `which is much slower but will give a nice plot in some cases.`, ``, `- Remaining arguments are interpreted as options which are`, `specified as equations of the form option = value. There are`, `currently 4 options -- grid, stepsize, iterations, and scheme.`, `The value for grid is a list, [m,n], where the current default is`, `[16,12]. The value for stepsize controls the step size in the`, `Runge/Kutta numerical scheme and determines how many points will`, `be plotted, the current default is (b-a)/20 and 20 points are`, `plotted. It is possible to decrease the stepsize without`, `increasing the number of points plotted (which slows down the`, `routine and leads to unreasonably large plot structures) by`, `setting the value for iterations. iterations = 10 decreases the`, `stepsize by a factor of 10, and every 10th point will be stored and`, `plotted. The option scheme=adaptive specifies that an adaptive`, `fourth-fifth Runge/Kutta numerical scheme will be used. This `, `should be used carefully because it will be much slower.`, ``, ``, `EXAMPLES: fieldplot((t,y)->t-y,-4..4,-3..3);`, `fieldplot((t,y)->y,0..1,0..3,{[0,1]});`, `fieldplot((t,y)->-t/y,-4..4,-4..4,{[0,1],[0,2],[0,3],[0,4]});`, `fieldplot((t,y)->-t/y,-4..4,-4..4,[0,1],scheme=adaptive);`, `fieldplot((t,y)->sin(t+y),-4..4,-4..4,grid=[16,14]);`, ``, `eq:=(t,y)->sin(t+y):`, `inits:={seq(seq([2*i,2*j],i=-2..2),j=-2..2)}:`, `fieldplot(eq,-4..4,-4..4,inits,grid=[0,0]);`, ``, `fieldplot((t,y)->-y^2+t,0..20,0..10,[0,1]);`, `fieldplot((t,y)->-y^2+t,0..20,0..10,[0,1],iterations=10);`, `fieldplot((t,y)->-y^2+t,0..20,0..10,[0,1],scheme=adaptive);`, ``, `eq4:=(x,y)->(1-x^2-y^2)/(y-x+2):`, `fieldplot(eq4,-2..2,-2..2,grid=[15,15]);`, ``, `eq4:=(x,y)->(1-x^2-y^2)/(y-x+2):`, `init4:=seq(seq([j,i],i=-2..2),j=-2..2):`, `fieldplot(eq4,-2..2,-2..2,{init4});`, ``, `eq4:=(x,y)->(1-x^2-y^2)/(y-x+2):`, `fieldplot(eq4,-2..2,-2..2,[1,1],scheme=adaptive);`, ``, `eq5:=(t,y)->y*tan(t);`, `init5:=seq(seq([i,j],i=-2..2),j=-2..2):`, `fieldplot(eq5,-6..6,-6..6,grid=[0,0],{init5},`, ` iterations=5,stepsize=0.1);`, ``, `eq6:=(t,y)->-y*tan(t);`, `init6:=seq(seq([2*i,2*j],i=-2..2),j=-2..2):`, `fieldplot(eq6,-6..6,-6..6,grid=[0,0],{init6},`, ` iterations=5,stepsize=.1);`, ` `, `see also: firsteuler, impeuler, rungekutta, rungekuttahf, phaseplot`): `help/text/phaseplot` := TEXT( `FUNCTION: phaseplot - draw phase space with flow lines for second`, ` order autonomous differential equation.`, ` `, `CALLING SEQUENCE:`, ` phaseplot([f1,f2], h, v)`, ` phaseplot([f1,f2], h, v,{inits},...)`, ``, `PARAMETERS: f1,f2 - two functions of three variables where`, `D(x)(t)=f1(t,x,y) and D(y)(t)=f2(t,x,y) where normally f1 and f2`, `do not depend on t. t=0 is substituted to draw the phase space`, `otherwise. Note it will still be interesting to draw flow`, `lines in the non-autonomous case, i.e. where f1 and f2 depend`, `explicity on time. h - horizontal range v - vertical range {inits}`, `- set of initial values for flow lines.`, ``, `SYNOPSIS: - A typical call to the phaseplot function is`, `phaseplot([f1,f2],a..b,c..d), where f1 and f2 are real functions of`, `three variables and a..b,c..d specifies the horizontal and`, `vertical real range. For each point, (xp,yp), of a grid, a scaled`, `vector with basepoint (xp,yp) and endpoint`, `(f1(0,xp,yp),f2(0,xp,yp)) is drawn.`, ``, `- Each initial point is specified as a list, [t0,x0,y0]. A`, `fourth order Runge/Kutta scheme is used to calculate the solution`, `to the system of equations, D(x)(t)=f1(t,x,y),`, `D(y)(t)=f2(t,x,y),x(t0)=x0,y(t0)=y0. The default is to use a step`, `size of 0.1 and do 50 iterations.`, ``, `- Remaining arguments are interpreted as options which are`, `specified as equations of the form option = value. There are`, `currently 4 options -- grid, stepsize, numsteps, and iterations.`, `The value for grid is a list, [m,n], where the current default is`, `[16,14]. The value for stepsize controls the step size in the`, `Runge/Kutta numerical scheme, the current default is 0.1. The`, `value for numsteps specifies how many steps forward and backward`, ` will be calculated and stored by the Runge/Kutta routine. The`, `current default is 50 forward and back. It is possible to`, `decrease the stepsize without increasing the number of points`, `plotted (which slows down the routine and leads to unreasonably`, `large plot structures) by setting the value for iterations.`, `iterations = 10, would decrease the stepsize by a factor of 10`, `and every 10th point would be stored and plotted.`, ` `, ``, `EXAMPLES: `, `phaseplot([(t,x,y)->-x,(t,x,y)->y],-4..4,-4..4);`, ``, `eq1:=(t,x,y)-> 5*x-x^2-x*y: `, `eq2:=(t,x,y)->x*y-2*y:`, `inits:={seq([0,5,i],i=1..9),seq([0,i,0.1],i=1..9)}:`, `phaseplot([eq1,eq2],0..10,0..10,inits,grid=[15,15]);`, ``, `readlib(Heaviside): `, `eq3:=(t,y,v)->v: `, `eq4:=(t,y,v)->-Heaviside(y):`, `phaseplot([eq3,eq4],-4..4,-4..4); `, ``, `readlib(Heaviside): `, `eq5:=(t,y,v)->v:`, `eq6:=(t,y,v)->-Heaviside(y-5)*v/2-(1-Heaviside(y-5))*v/10+1:`, `phaseplot([eq5,eq6],0..10,0..5,[0,0,0],iterations=10,`, ` numsteps=100,grid=[15,15]);`, ``, `eq7:=(t,x,y)-> -x+x*y: `, `eq8:=(t,x,y)->-y+2*x*y:`, `inits2:={seq(seq([0,i/2,j/2],i=0..2),j=0..2)}:`, `phaseplot([eq7,eq8],-.5..1.5,-.5..1.5,inits2,grid=[15,15]);`, ``, `pend3:=(t,y,v) -> v;`, `pend4:=(t,y,v) -> -sin(y) -0.2*v;`, `pendinits:=seq([0,j/2],j=2..8):`, `phaseplot([pend3,pend4],-10..10,-4..4,{pendinits},`, ` numsteps=200);`, ``, `ropedv:= proc(t,v,y) if v > 0 then (9800 -3*980*y - 3*v^2)/(3*y)`, ` else (9800 -3*980*y)/(3*y) fi end;`, `ropedy:= (t,v,y) -> v;`, `phaseplot([ropedv,ropedy],-59..59,0.1..5);`, ``, `eq11:=(t,y,v)->v^2-y^2;`, `eq12:=(t,y,v)->y-2*v;`, `eq1init:=seq(seq([i,j],j=-2..2),i=-2..2):`, `phaseplot([eq11,eq12],-5..5,-5..5,grid=[12,12],{eq1init});`, ``, `eq21:=(t,x,y)->cos(x+y);`, `eq22:=(t,x,y)->-sin(y);`, `phaseplot([eq21,eq22],-6..6,-6..6,grid=[25,25]);`, ``, `see also: firsteuler, impeuler, rungekutta, rungekuttahf, fieldplot` ): `help/text/rungekutta` := TEXT( `FUNCTION: rungekutta - approximate a solution to a system of first`, `order ordinary differential equations using a fourth order Runge/Kutta`, `numerical method. Solution is returned as an array of lists.`, `rungekuttahf is an alternative form of this procedure which uses`, `hardware floating point operations and is not affected by the Digits`, `settings. rungekuttahf is limited to double precision on most systems`, `and will be about as accurate as rungekutta using a Digits setting `, `of 15. It is however, MUCH faster.`, ``, `CALLING SEQUENCE:`, ` rungekutta(fcnlist,initvalue,stepsize,numsteps)`, ` rungekutta(fcnlist,initvalue,stepsize,numsteps,iterations)`, ``, `PARAMETERS:`, ` fcnlist - list of functions, [f1,f2,...,fn], of n+1 real variables. `, ` The differential equation to be solved is of the form:`, ` dx1/dt = f1(t,x1,x2,...,xn)`, ` dx2/dt = f2(t,x1,x2,...,xn)`, ` . . .`, ` . . .`, ` . . .`, ` dxn/dt = fn(t,x1,x2,...,xn)`, ` initvalue - list of n+1 numbers, [i0,i1,i2,...,in], which represents `, ` the initial value of the differential equation. `, ` stepsize - step size used by the numerical method.`, ` numsteps - number of steps taken by numerical method.`, ` iterations - optional value which decreases the step size without `, ` storing more points. New step size will be `, ` stepsize/iterations.`, ``, `SYNOPSIS: - A typical call to the rungekutta function is`, `rungekutta([f1,f2,...,fn],[i0,i1,i2,...,in],stepsize,numsteps), where`, `f1,f2,...,fn are real valued functions of n+1 variables, stepsize is`, `the step size used by the numerical method, and numsteps specifies the`, `number of steps to be taken by the numerical method. An array of (n+1)`, `lists of (numsteps+1) numbers is returned, where each list represents`, `a point calculated by the numerical method.`, ``, `- An optional fifth argument is used to reduce the step size without `, `storing more points. This is especially useful when graphing `, `solutions which require a small step size and a large plot structure `, `is undesirable. Also useful when only the final value is wanted.`, ``, `- makelist is a procedure in the ODE file to extract a list of values`, `from an array of lists. This is useful for making two dimensional`, `plots of solutions calculated by rungekutta. makelist(A,m,n) will `, `make a list of values from the mth and nth components from each `, `entry of A. m defaults to 1 and n defaults to 2 if no values are `, `entered.`, ``, `EXAMPLES: `, `rkpts1:=rungekutta((t,y)->t-y,[0,1],.1,10);`, `plot(makelist(rkpts1));`, ``, `eq1:=(t,x,y)-> y;`, `eq2:=(t,x,y)->-x;`, `rkpts2:=rungekutta([eq1,eq2],[0,1,1],.1,60);`, `plot(makelist(rkpts2,2,3));`, ``, `dp1:=(t,x,xp,y,yp)->xp;`, `dp2:=(t,x,xp,y,yp)->-2*sin(x)+sin(y);`, `dp3:=(t,x,xp,y,yp)->yp;`, `dp4:=(t,x,xp,y,yp)->-2*sin(y)+2*sin(x);`, `dpinit:=[0,0,0,0,2];`, `dppts:=rungekutta([dp1,dp2,dp3,dp4],dpinit,1,10,1):`, `dppts2:=rungekutta([dp1,dp2,dp3,dp4],dpinit,1,10,2):`, `dppts3:=rungekutta([dp1,dp2,dp3,dp4],dpinit,1,10,4):`, `plot({makelist(dppts,2,4),makelist(dppts2,2,4),`, ` makelist(dppts3,2,4)});`, ``, `see also: firsteuler, impeuler, rungekuttahf, fieldplot, phaseplot`): `help/text/rungekuttahf` := TEXT( `FUNCTION: rungekuttahf - approximate a solution to a system of up to`, `four first order ordinary differential equations using a fourth order`, `Runge/Kutta numerical method. Solution is returned as an array of`, `lists. rungekutta is an alternative form of this procedure which uses`, `Maple's arbritrary precision floating point operations and is thus`, `affected by the Digits settings. rungekutta is useful for studying the`, `effect of the round off error that occurs with numerical methods or`, `with differential equations for which rungekuttahf fails. rungekuttahf`, `is limited to double precision on most systems and will be about as`, `accurate as rungekutta using a Digits setting of 15. rungekuttahf is`, `however, MUCH faster.`, ``, `CALLING SEQUENCE:`, ` rungekuttahf(fcnlist,initvalue,stepsize,numsteps)`, ` rungekuttahf(fcnlist,initvalue,stepsize,numsteps,iterations)`, ``, `PARAMETERS:`, ` fcnlist - list of functions, [f1,f2,...,fn], of n+1 real variables. `, ` The differential equation to be solved is of the form:`, ` dx1/dt = f1(t,x1,x2,...,xn)`, ` dx2/dt = f2(t,x1,x2,...,xn)`, ` . . .`, ` . . .`, ` . . .`, ` dxn/dt = fn(t,x1,x2,...,xn)`, ` initvalue - list of n+1 numbers, [i0,i1,i2,...,in], which represents `, ` the initial value of the differential equation. `, ` stepsize - step size used by the numerical method.`, ` numsteps - number of steps taken by numerical method.`, ` iterations - optional value which decreases the step size without `, ` storing more points. New step size will be `, ` stepsize/iterations.`, ``, `SYNOPSIS: - A typical call to the rungekuttahf function is`, `rungekuttahf([f1,f2,...,fn],[i0,i1,i2,...,in],stepsize,numsteps), where`, `f1,f2,...,fn are real valued functions of n+1 variables, stepsize is`, `the step size used by the numerical method, and numsteps specifies the`, `number of steps to be taken by the numerical method. An array of (n+1)`, `lists of (numsteps+1) numbers is returned, where each list represents`, `a point calculated by the numerical method.`, ``, `- An optional fifth argument is used to reduce the step size without `, `storing more points. This is especially useful when graphing `, `solutions which require a small step size and a large plot structure `, `is undesirable. Also useful when only the final value is wanted.`, ``, `- makelist is a procedure in the ODE file to extract a list of values`, `from an array of lists. This is useful for making two dimensional`, `plots of solutions calculated by rungekuttahf. makelist(A,m,n) will `, `make a list of values from the mth and nth components from each `, `entry of A. m defaults to 1 and n defaults to 2 if no values are `, `entered.`, ``, `EXAMPLES: `, `rkpts1:=rungekuttahf((t,y)->t-y,[0,1],.1,10);`, `plot(makelist(rkpts1));`, ``, `eq1:=(t,x,y)-> y;`, `eq2:=(t,x,y)->-x;`, `rkpts2:=rungekuttahf([eq1,eq2],[0,1,1],.1,60);`, `plot(makelist(rkpts2,2,3));`, ``, `dp1:=(t,x,xp,y,yp)->xp;`, `dp2:=(t,x,xp,y,yp)->-2*sin(x)+sin(y);`, `dp3:=(t,x,xp,y,yp)->yp;`, `dp4:=(t,x,xp,y,yp)->-2*sin(y)+2*sin(x);`, `dpinit:=[0,0,0,0,2];`, `dppts:=rungekuttahf([dp1,dp2,dp3,dp4],dpinit,.1,100):`, `plot({makelist(dppts,2,4)});`, ``, `eq3:=(t,x,y)->x*sin(t)+5*y;`, `eq4:=(t,x,y)->-10*x+y*sin(t);`, `rkpts2:=rungekuttahf([eq3,eq4],[0,2,0],.1,100,5);`, `plot(makelist(rkpts2,2,3));`, `plot({makelist(rkpts2,1,2),makelist(rkpts2,1,3)});`, `#for the Maple expert and the brave alike, construct a PLOT3D structure`, `plot1:=PLOT3D(CURVELIST([seq([rkpts2[i][1],rkpts2[i][2],rkpts2[i][3]],`, ` i=0..100)]),AXES(BOXED)):`, `plot1;`, ` `, ``, `see also: firsteuler, impeuler, rungekutta, fieldplot, phaseplot`): `help/text/firsteuler` := TEXT( `FUNCTION: firsteuler - approximate a solution to a first order`, `ordinary differential equations using Euler's method.`, `Solution is returned as an array of lists. `, ``, `CALLING SEQUENCE:`, ` firsteuler(f,i,h,n)`, ``, `PARAMETERS:`, ` f - function of 2 real variables. The differential equation to be `, ` solved is dy/dt = f(t,y) `, ` i - list of 2 numbers. `, ` h - step size used by the numerical method.`, ` n - number of steps taken by numerical method.`, ``, `SYNOPSIS: - A typical call to the firsteuler function is`, `firsteuler(f,i,h,n) where f is a real valued functions of 2 variables,`, `h is the step size used by the numerical method, and n specifies the`, `number of steps to be taken by the numerical method. An array of (n+1)`, `lists of 2 numbers is returned, where each list represents a point`, `calculated by the numerical method.`, ``, `- makelist is a procedure in the ODE file to extract a list of values`, `from an array of lists. This is useful for making a plot of the`, `approximate solution calculated by firsteuler. makelist(A) will make a`, `list of values from the entries of A. `, ``, `EXAMPLES: `, `eulerpts:=firsteuler((t,y)->t-y,[0,1],.1,10);`, `plot(makelist(eulerpts));`, ``, `eq:=(t,y)->y;`, `pts1:=firsteuler(eq,[0,1],.5,2);`, `pts2:=firsteuler(eq,[0,1],.25,4);`, `pts3:=firsteuler(eq,[0,1],.1,10);`, `plot({makelist(pts1),makelist(pts2),makelist(pts3),`, ` exp(t)},0..1,style=LINE);`, ` `, `eq2:=(t,y)->-t/y;`, `eq2pts:=firsteuler(eq2,[0,4],.2,25):`, `plot(makelist(eq2pts));`, ``, `eq3:=(t,y)->5*y-6*exp(-t);`, `eq3pts:=firsteuler(eq3,[0,1],.1,10);`, `eq3sol:=dsolve({diff(y(t),t)=eq3(t,y),y(0)=1},y(t));`, `plot({makelist(eq3pts),rhs(eq3sol)},0..1);`, ``, `eq4:=(t,y)->25*y*(1-y);`, `eq4pts:=firsteuler(eq4,[0,1.3],.1,10):`, `plot(makelist(eq4pts));`, ``, `eq4:=(t,y)->25*y*(1-y);`, `eq4pts:=firsteuler(eq4,[0,1.3],.1,10):`, `eq4pts2:=firsteuler(eq4,[0,1.3],.05,20):`, `eq4pts3:=firsteuler(eq4,[0,1.3],.025,40):`, `plot({makelist(eq4pts),makelist(eq4pts2),makelist(eq4pts3)});`, ``, ``, `see also: impeuler, rungekutta, rungekuttahf, fieldplot, phaseplot`): `help/text/impeuler` := TEXT( `FUNCTION: impeuler - approximate a solution to a first order ordinary`, `differential equations using the improved Euler method. Solution is`, `returned as an array of lists. `, ``, `CALLING SEQUENCE:`, ` impeuler(f,i,h,n)`, ``, `PARAMETERS:`, ` f - function of 2 real variables. The differential equation to be `, ` solved is dy/dt = f(t,y) `, ` i - list of 2 numbers. `, ` h - step size used by the numerical method.`, ` n - number of steps taken by numerical method.`, ``, `SYNOPSIS: - A typical call to the impeuler function is`, `impeuler(f,i,h,n) where f is a real valued functions of 2 variables,`, `h is the step size used by the numerical method, and n specifies the`, `number of steps to be taken by the numerical method. An array of (n+1)`, `lists of 2 numbers is returned, where each list represents a point`, `calculated by the numerical method.`, ``, `- makelist is a procedure in the ODE file to extract a list of values`, `from an array of lists. This is useful for making a plot of the`, `approximate solution calculated by impeuler. makelist(A) will make a`, `list of values from the entries of A. `, ``, `EXAMPLES: `, `impts:=impeuler((t,y)->t-y,[0,1],.1,10);`, `plot(makelist(impts));`, ``, `eq:=(t,y)->y;`, `impts1:=impeuler(eq,[0,1],.5,2);`, `impts2:=impeuler(eq,[0,1],.25,4);`, `impts3:=impeuler(eq,[0,1],.1,10);`, `plot({makelist(impts1),makelist(impts2),makelist(impts3),`, ` exp(t)},0..1,style=LINE);`, ` `, `eq2:=(t,y)->-t/y;`, `eq2impts:=impeuler(eq2,[0,4],.2,25):`, `plot(makelist(eq2impts));`, ``, `eq3:=(t,y)->5*y-6*exp(-t);`, `eq3impts:=impeuler(eq3,[0,1],.1,10);`, `eq3sol:=dsolve({diff(y(t),t)=eq3(t,y),y(0)=1},y(t));`, `plot({makelist(eq3impts),rhs(eq3sol)},0..1);`, ``, `eq4:=(t,y)->25*y*(1-y);`, `eq4impts:=impeuler(eq4,[0,1.3],.1,10):`, `plot(makelist(eq4impts));`, ``, `eq4:=(t,y)->25*y*(1-y);`, `eq4impts:=impeuler(eq4,[0,1.3],.1,10):`, `eq4impts2:=impeuler(eq4,[0,1.3],.05,20):`, `eq4impts3:=impeuler(eq4,[0,1.3],.025,40):`, `eq4implot:=plot({makelist(eq4impts),makelist(eq4impts2),`, ` makelist(eq4impts3)}):`, `eq4field:=fieldplot(eq4,0..1,0.8..1.5):`, `with(plots):`, `display({eq4implot,eq4field});`, ``, ``, `see also: firsteuler, rungekutta, rungekuttahf, fieldplot, phaseplot`): ################################################# firsteuler := proc(f,initlist,h0,n,iterations) local i,xk,yk,pts,h,npts; options `Copyright 1991 by Daniel Schwalbe`; pts:=array(0..n); if nargs>4 then npts:=iterations; else npts:=1 fi; h:=evalf(h0/npts); xk:=evalf(initlist[1]); yk:=evalf(initlist[2]); pts[0] := [xk,yk]; for i to n while type(yk,numeric) do to npts while type(yk,numeric) do yk:=yk + h*f(xk,yk); xk:=xk+h od; pts[i]:=[xk,yk] od: if i <= n then ERROR(`undefined procedure`) else op(pts) fi: end: ################################################## impeuler:=proc(f,initlist,h0,n,iterations) local h,xk,yk,pts,j,npts: options `Copyright 1991 by Daniel Schwalbe`; pts:=array(0..n); if nargs>4 then npts:=iterations else npts:=1 fi; h:=evalf(h0/npts): xk:=evalf(initlist[1]); yk:=evalf(initlist[2]); pts[0] := [xk,yk]; for j to n while type(yk,numeric) do to npts while type(yk,numeric) do yk:=yk+h/2*(f(xk,yk)+f(xk+h,yk+h*f(xk,yk))): xk:=xk+h od: pts[j]:=[xk,yk] od: if j <= n then ERROR(`undefined procedure`) else op(pts) fi: end: ################################################## rungekutta:=proc(f,initlist,h0,n,iterations) local k1,k2,k3,k4,fl,h,pts,j,i,k,m,npts,temp: options `Copyright 1991 by Daniel Schwalbe`; m:=nops([op(initlist)])-1; fl:=f; fl:=[op(fl)]; pts:=array(0..n); if nargs > 4 then npts:=iterations else npts:=1 fi; h:=evalf(h0/npts): pts[0] := map(evalf,[op(initlist)]); temp:=pts[0]; for j from 0 to n-1 while not member(false,map(type,(temp,numeric))) do to npts while not member(false,map(type,(temp,numeric))) do for i to m do k1[i]:= evalf(fl[i](op(temp))) od: for i to m do k2[i]:= evalf(fl[i](temp[1]+h/2, seq(temp[k+1] + h*k1[k]/2,k=1..m))) od: for i to m do k3[i]:=evalf(fl[i](temp[1]+h/2, seq(temp[k+1] + h*k2[k]/2,k=1..m))) od: for i to m do k4[i]:=evalf(fl[i](pts[j][1]+h, seq(temp[k+1] + h*k3[k],k=1..m))) od: temp:=[temp[1]+h, seq(temp[i+1]+ h*(k1[i]+2*k2[i]+2*k3[i]+k4[i])/6,i=1..m)] od: pts[j+1]:= temp od: if j<=n-1 then ERROR(`undefined procedure`) else op(pts) fi: end: ################################################## rkhf:=proc(Akf,B,h,n,npts) local j,i,AA,k: options `Copyright 1991 by Daniel Schwalbe`; AA:=array(0..npts,1..2); for i to 2 do AA[0,i]:=B[0,i] od; for i to n do for j to npts do AA[j,1]:=AA[j-1,1]+h; AA[j,2]:=Akf(AA,j-1) od; for k to 2 do B[i,k]:=AA[npts,k] od; for k to 2 do AA[0,k]:=AA[npts,k] od; od; end: rkhf2:=proc(Akf,Akf2,B,h,n,npts) local j,i,AA,k: options `Copyright 1991 by Daniel Schwalbe`; AA:=array(0..npts,1..3); for i to 3 do AA[0,i]:=B[0,i] od; for i to n do for j to npts do AA[j,1]:=AA[j-1,1]+h; AA[j,2]:=Akf(AA,j-1); AA[j,3]:=Akf2(AA,j-1) od: for k to 3 do B[i,k]:=AA[npts,k] od; for k to 3 do AA[0,k]:=AA[npts,k] od; od; end: rkhf3:=proc(Akf,Akf2,Akf3,B,h,n,npts) local j,i,AA,k: options `Copyright 1991 by Daniel Schwalbe`; AA:=array(0..npts,1..4); for i to 4 do AA[0,i]:=B[0,i] od; for i to n do for j to npts do AA[j,1]:=AA[j-1,1]+h; AA[j,2]:=Akf(AA,j-1); AA[j,3]:=Akf2(AA,j-1): AA[j,4]:=Akf3(AA,j-1) od: for k to 4 do B[i,k]:=AA[npts,k] od; for k to 4 do AA[0,k]:=AA[npts,k] od; od; end: rkhf4:=proc(Akf,Akf2,Akf3,Akf4,B,h,n,npts) local j,i,AA,k: options `Copyright 1991 by Daniel Schwalbe`; AA:=array(0..npts,1..5); for i to 5 do AA[0,i]:=B[0,i] od; for i to n do for j to npts do AA[j,1]:=AA[j-1,1]+h; AA[j,2]:=Akf(AA,j-1); AA[j,3]:=Akf2(AA,j-1): AA[j,4]:=Akf3(AA,j-1): AA[j,5]:=Akf4(AA,j-1) od: for k to 5 do B[i,k]:=AA[npts,k] od; for k to 5 do AA[0,k]:=AA[npts,k] od; od; end: rungekuttahf:=proc(f,init,h0,n,iterations) local h,i,m,A,B,fl,k,k1,k2,k3,k4,precision, npts,Ak,Akf,Akf2,Akf3,Akf4,j; options `Copyright 1991 by Daniel Schwalbe`; precision := max(Digits,trunc(evalhf(Digits))); m:=nops([op(init)])-1; A:=array(0..n,1..m+1); if nargs > 4 then npts:=iterations else npts:=1 fi; h:=evalf(h0/npts,precision): for i to m+1 do A[0,i]:=evalf(init[i],precision) od; if type(f,list) then fl:=f else fl:=[f] fi; for i to m do k1[i]:=fl[i](seq(A[j,k],k=1..m+1)) od: for i to m do k2[i]:=fl[i](A[j,1]+h/2, seq(A[j,k+1] + h*k1[k]/2,k=1..m)) od: for i to m do k3[i]:=fl[i](A[j,1]+h/2, seq(A[j,k+1] + h*k2[k]/2,k=1..m)) od: for i to m do k4[i]:=fl[i](A[j,1]+h, seq(A[j,k+1] + h*k3[k],k=1..m)) od: for i to m do Ak[i]:= A[j,i+1]+h*(k1[i]+2*k2[i]+2*k3[i]+k4[i])/6 od: if m=1 then Akf:=unapply(Ak[1],A,j); evalhf(rkhf(Akf,var(A),h,n,npts)) elif m=2 then Akf:=unapply(Ak[1],A,j); Akf2:=unapply(Ak[2],A,j); evalhf(rkhf2(Akf,Akf2,var(A),h,n,npts)) elif m=3 then Akf:=unapply(Ak[1],A,j); Akf2:=unapply(Ak[2],A,j); Akf3:=unapply(Ak[3],A,j); evalhf(rkhf3(Akf,Akf2,Akf3,var(A),h,n,npts)) elif m=4 then Akf:=unapply(Ak[1],A,j); Akf2:=unapply(Ak[2],A,j); Akf3:=unapply(Ak[3],A,j); Akf4:=unapply(Ak[4],A,j); evalhf(rkhf4(Akf,Akf2,Akf3,Akf4,var(A),h,n,npts)) else print(`rungekuttahf is not designed for more than 4 equations.`); print(`Try rungekutta as an alternate.`) fi: B:=array(0..n); for i from 0 to n do B[i]:=[seq(A[i,k],k=1..m+1)] od; op(B) end: ################################################## makelist:=proc(A,c10,c20) local c1,c2,ri,rf,i; options `Copyright 1991 by Daniel Schwalbe`; if nargs = 1 then ri:=op(1,op(2,eval(A))); rf:=op(2,op(2,eval(A))); c1:=1; c2:=2; else c1:=c10; c2:=c20; ri:=op(1,op(2,eval(A))); rf:=op(2,op(2,eval(A))) fi; map(op,[seq([[A[i][c1],A[i][c2]]],i=ri..rf)]); end: ################################################# intbox:=proc(p,q,r,s,pta,pt) local xnew,ynew; options `Copyright 1991 by Daniel Schwalbe`; if pt[1]

q then xnew:=q; ynew:=(pt[2]-pta[2])/(pt[1]-pta[1])*(xnew-pta[1])+pta[2]; pt[1]:=xnew; pt[2]:=ynew fi; if pt[2]s then ynew:=s; xnew:=(ynew-pta[2])*(pt[1]-pta[1])/(pt[2]-pta[2])+pta[1]; pt[1]:=xnew; pt[2]:=ynew fi; end: ################################################## rkstep:=proc(f,hp,endpt,npts,pt) local h,t,y,k1,k2,k3,k4; options `Copyright 1991 by Daniel Schwalbe`; t:=pt[1]:y:=pt[2]: h:=hp: to npts while ((h>0 and tendpt)) do if (h>0 and t+h>endpt) or (h<0 and t+h=r do A[2*j+1]:=pt[2]: A[2*j]:=pt[1]; CC[9]:=2*j+1: j:=j+1; rkstep(f,h,q,npts,pt): od: if j>0 then pta[1]:=A[2*j-2]; pta[2]:=A[2*j-1]; intbox(p,q,r,s,pta,pt); A[2*j]:=pt[1]; A[2*j+1]:=pt[2]; CC[9]:=2*j+1 fi; j:=0;CC[8]:=0; pt[1]:=CC[1]; pt[2]:=CC[2]; while pt[1] > p and pt[2]>=r and pt[2]<=s do A[2*j+1]:=pt[2]: A[2*j]:=pt[1]; CC[8]:=2*j: j:=j-1; rkstep(f,-h,p,npts,pt): od: if j<0 then pta[1]:=A[2*j+2]; pta[2]:=A[2*j+3]; intbox(p,q,r,s,pta,pt); A[2*j]:=pt[1]; A[2*j+1]:=pt[2]; CC[8]:=2*j fi; end: adaptiverk:=proc(f,CC) local F,j,tt,Fval,yy,A,p,q,h,t,y; options `Copyright 1991 by Daniel Schwalbe`; p:=CC[4]: q:=CC[5]: tt:=CC[1]:yy:=CC[2]:h:=CC[3]: A[0]:=tt,yy: F:=dsolve({D(y)(t)=f(t,y(t)),y(tt)=yy},y(t),numeric); j:=0; Fval:=traperror(F(tt+h)); while Fval<>lasterror and tt+h<=q do tt:=tt+h; j:=j+1; A[j]:=Fval: Fval:=traperror(F(tt+h)) od; CC[9]:=j: j:=0;tt:=CC[1]: Fval:=traperror(F(tt-h)); while Fval<>lasterror and tt-h>=p do tt:=tt-h; j:=j-1; A[j]:=Fval: Fval:=traperror(F(tt-h)) od; CC[8]:=j: op(A) end: computelineshf:=proc(f,Ax,Ay,CC) local m,n,dx,dy,k,slp,xinc,yinc,a,b,p,q,r,s,i,j,ldx,ldy; m:=CC[1]; n:=CC[2]; p:=CC[4]; q:=CC[5]; r:=CC[6]; s:=CC[7]; dx:=(q-p)/m; dy:=(s-r)/n; ldx:=3*dx/8; ldy:=3*dy/8; a:=p; for i from 0 to m do b:=r; for j from 0 to n do k:=2*j; slp:=f(a,b); if slp > dy/dx or -slp > dy/dx then xinc:=ldy/slp; yinc:=ldy else xinc:=ldx; yinc:=ldx*slp fi; Ax[i,k]:=a-xinc; Ay[i,k]:=b-yinc; Ax[i,k+1]:=a+xinc; Ay[i,k+1]:=b+yinc; b:=b+dy od; a:=a+dx od; end: computelines:=proc(f,Ax,Ay,p,q,r,s,m,n) local dx,dy,k,ldx,slp,ldy,xinc,yinc,a,b,i,j; dx:=(q-p)/m; dy:=(s-r)/n; ldx:=3*dx/8; ldy:=3*dy/8; a:=p; for i from 0 to m do b:=r; for j from 0 to n do k:=2*j; slp:=traperror(evalf(f(a,b))); if lasterror=slp then xinc:=0; yinc:=0; elif not type(slp,'numeric') then xinc:=0; yinc:=0; elif slp > dy/dx or -slp > dy/dx then xinc:=ldy/slp; yinc:=ldy; else xinc:=ldx; yinc:=ldx*slp fi; Ax[i,k]:=a-xinc; Ay[i,k]:=b-yinc; Ax[i,k+1]:=a+xinc; Ay[i,k+1]:=b+yinc; b:=b+dy od; a:=a+dx od; end: macro( INITCOND1 = {[constant,constant]} ); fieldplot:=proc(F,H,V) local npts,init,flow,A,Ax,Ay,i,j,nsteps,f,numstyle, CC,h,m,n,p,q,r,s,cargs,flows,precision,step, fl, flines,carg,initlist,Aer: options `Copyright 1991 by Daniel Schwalbe`; precision := max(Digits,trunc(evalhf(Digits))); npts:=1; m:=16; n:=12; nsteps:=20: if not type(F,{procedure,numeric}) then ERROR(`1st argument (functions) should be a procedure`) fi; if not type([H,V],[constant..constant,constant..constant]) then ERROR(`2nd and 3rd arguments (domains) must be real ranges`) fi; step:= proc(x) if x<0 then 0 else 1 fi end; if type(F,numeric) then f:=subs(_body=F,proc() _body end) else f:=subs(Heaviside=step,eval(F)) fi: p:=evalf(lhs(H),precision); q:=evalf(rhs(H),precision); r:=evalf(lhs(V),precision); s:=evalf(rhs(V),precision); h:=evalf((q-p)/nsteps,precision): initlist:={}; cargs:={args[4..nargs]}; for carg in cargs do if type(carg,'set') then if type(carg,'set(INITCOND1)') then initlist := carg else ERROR(`initial conditions must be a set of points of type`,INITCOND1) fi; elif type(carg,'list') then # A single initial condition if type(carg,INITCOND1) then initlist:={carg} else ERROR(`initial condition must be of type`,INITCOND1) fi; elif not type(carg,equation) then ERROR(`optional arguments must be equations`) elif op(1,carg)='stepsize' then h:=evalf(op(2,carg),precision); nsteps:=trunc((q-p)/h + 1); if not type(h,numeric) then ERROR(`stepsize must be a numerical value`) fi; elif op(1,carg)='numsteps' then nsteps:=op(2,carg); if not type(nsteps,posint) then ERROR(`numsteps must be a positive integer`) fi; elif op(1,carg)='iterations' then npts:=op(2,carg); if not type(npts,posint) then ERROR(`number of iterations must be a positive integer`) fi; elif op(1,carg)='scheme' then numstyle := 'adaptive' elif op(1,carg)='grid' then if not type(rhs(carg),[integer,integer]) then ERROR(`grid must be a list of two integers`) fi; m:=op(2,carg)[1]; n:=op(2,carg)[2] else ERROR(`unknown option`,lhs(carg)) fi od; CC:=array(1..9); CC[1]:=m: CC[2]:=n: CC[3]:=h: CC[4]:=p: CC[5]:=q: CC[6]:=r: CC[7]:=s: if m>0 and n>0 then Ax:=array(0..m,0..2*n+1); Ay:=array(0..m,0..2*n+1); fl:=traperror(evalhf(computelineshf( f,var(Ax),var(Ay),var(CC)))): if lasterror=fl then #print(`computelineshf failed`); computelines(f,Ax,Ay,p,q,r,s,m,n) fi; flines:=seq(seq(CURVE(BLACK,LINE, [Ax[i,2*j],Ay[i,2*j],Ax[i,2*j+1], Ay[i,2*j+1]]),i=0..m),j=0..n): else flines:=[0,0] fi: i:=1: for init in initlist do CC[1]:=evalf(init[1],precision); CC[2]:=evalf(init[2],precision); if numstyle<>'adaptive' then A:=array(2*(-nsteps-3)..2*(nsteps)+3): Aer:=traperror(evalhf(fieldrk(var(A),f,var(CC),npts))); if Aer=lasterror then #print(`evalhf(fieldrk) failed`); Aer:= traperror(fieldrk(A,f,CC,npts)) fi; if Aer<>lasterror then flow[i]:=seq(A[j],j=round(CC[8])..round(CC[9])): i:=i+1 fi; else A:=adaptiverk(f,CC): flow[i]:=seq(A[j],j=round(CC[8])..round(CC[9])); i:=i+1 fi od: flows:=seq(CURVE(RED,LINE,[flow[j]]),j=1..i-1); m:=max(m,16);n:=max(n,12); PLOT(AXIS(HORIZONTAL,SCALE(4,LINEAR,NUMBER)), AXIS(VERTICAL,SCALE(4,LINEAR,NUMBER)), RANGE(HORIZONTAL,p-(q-p)/m..q+(q-p)/m), RANGE(VERTICAL,r-(s-r)/n..s+(s-r)/n), flines,flows): end: ################################################## rk2step:=proc(f1,f2,h,p,q,r,s,npts,pt) local t,x,y,k11,k12,k13,k14,k21,k22,k23,k24; options `Copyright 1991 by Daniel Schwalbe`; t:=pt[1]:x:=pt[2]:y:=pt[3]: to npts while x<=q and x>=p and y<=s and y>=r do k11:=f1(t,x,y): k21:=f2(t,x,y): k12:=f1(t+h/2,x+h*k11/2,y + h*k21/2): k22:=f2(t+h/2,x+h*k11/2,y + h*k21/2): k13:=f1(t+h/2,x+h*k12/2,y + h*k22/2): k23:=f2(t+h/2,x+h*k12/2,y + h*k22/2): k14:=f1(t+h,x+h*k13,y + h*k23): k24:=f2(t+h,x+h*k13,y + h*k23): x:=x+h*(k11+2*k12+2*k13+k14)/6; y:=y+h*(k21+2*k22+2*k23+k24)/6; t:=t+h; od; pt[1]:=t: pt[2]:=x:pt[3]:=y: end: phaserk:=proc(A,f1,f2,CC,nsteps,npts) local pta,ptb,pt,p,q,r,s,j,h; options `Copyright 1991 by Daniel Schwalbe`; pt:=array(1..3): pta:=array(1..2):ptb:=array(1..2): h:=CC[1]/npts: pt[1]:=CC[2]: pt[2]:=CC[3]: pt[3]:=CC[4]: A[0]:=pt[2] :A[1]:=pt[3]: p:=CC[5]:q:=CC[6]:r:=CC[7]:s:=CC[8]: j:=0: while j=p and pt[3]<=s and pt[3]>=r do j:=j+1: rk2step(f1,f2,h,p,q,r,s,npts,pt): A[2*j]:=pt[2]: A[2*j+1]:=pt[3]: od: if pt[2]>q or pt[2]

s or pt[3]=p and pt[3]<=s and pt[3]>=r do j:=j-1: rk2step(f1,f2,-h,p,q,r,s,npts,pt): A[2*j]:=pt[2]: A[2*j+1]:=pt[3]: od: if pt[2]>q or pt[2]

s or pt[3] n then blocks:= m else blocks:= n fi; maxlen:=0; a:=p; for j from 0 to m do b:=r; for i from 0 to n do x:=fa(0,a,b); y:=fb(0,a,b); len:=((blocks/(q-p)*x)^2+(blocks/(s-r)*y)^2)^(1/2); if len > maxlen then maxlen:=len fi; b:=b+dy od; a:=a+dx od; a:=p; if blocks < 20 then blocks:=20 fi; for j from 0 to m do b:=r; for i from 0 to n do k:=6*i; x:=fa(0,a,b)/maxlen; y:=fb(0,a,b)/maxlen; if x=0 and y=0 then scale:= 0 else scale:=(((blocks*x/(q-p))^2+(blocks*y/(s-r))^2))^(-1/2) fi; wx:= scale*x; wy:=scale*y; #zx and zy are components of perpendicular vector zx:= -scale*y/(s-r)*(q-p); zy:= scale*x/(q-p)*(s-r) ; Ax[j,k]:=a; Ay[j,k]:=b; ap:=a+2*x/3+wx/3; bp:=b+2*y/3+wy/3; Ax[j,k+1]:=ap; Ay[j,k+1]:=bp; Ax[j,k+2]:=ap-wx/6+zx/6; Ay[j,k+2]:=bp-wy/6+zy/6; Ax[j,k+3]:=ap;#+wx/6; Ay[j,k+3]:=bp;#+wy/6; Ax[j,k+4]:=ap-wx/6-zx/6; Ay[j,k+4]:=bp-wy/6-zy/6; Ax[j,k+5]:=ap; Ay[j,k+5]:=bp; b:=b+dy od; a:=a+dx od end: #Only difference in computevects from computevectshf is trapping #for errors when evaluating input procedures. computevects:=proc(Ax,Ay,fa,fb,p,q,r,s,m,n) local i,j,k,maxlen,len,a,b,x,y,ap,bp, scale,wx,wy,zx,zy,blocks,dx,dy; options `Copyright 1991 by Daniel Schwalbe`; dx:=(q-p)/m: dy:=(s-r)/n: if m > n then blocks:= m else blocks:= n fi; maxlen:=0; a:=p; for j from 0 to m do b:=r; for i from 0 to n do x:= traperror(fa(0,a,b)); if x = lasterror then y:=0; x:=0 else y:=traperror(fb(0,a,b)) fi; if y = lasterror then y:=0; x:=0 fi; if not type([x,y],[numeric,numeric]) then x:=0; y:=0 fi; len:=((blocks/(q-p)*x)^2+(blocks/(s-r)*y)^2)^(1/2); if len > maxlen then maxlen:=len fi; b:=b+dy od; a:=a+dx od; a:=p; if blocks < 20 then blocks:=20 fi; for j from 0 to m do b:=r; for i from 0 to n do k:=6*i; x:= traperror(fa(0,a,b)/maxlen); if x = lasterror then y:=0; x:=0 else y:=traperror(fb(0,a,b)/maxlen) fi; if y = lasterror then y:=0; x:=0 fi; if not type([x,y],[numeric,numeric]) then x:=0; y:=0 fi; if x=0 and y=0 then scale:= 0 else scale:=(((blocks*x/(q-p))^2+(blocks*y/(s-r))^2))^(-1/2) fi; wx:= scale*x; wy:=scale*y; #zx and zy are components of perpendicular vector zx:= -scale*y/(s-r)*(q-p); zy:= scale*x/(q-p)*(s-r) ; Ax[j,k]:=a; Ay[j,k]:=b; ap:=a+2*x/3+wx/3; bp:=b+2*y/3+wy/3; Ax[j,k+1]:=ap; Ay[j,k+1]:=bp; Ax[j,k+2]:=ap-wx/6+zx/6; Ay[j,k+2]:=bp-wy/6+zy/6; Ax[j,k+3]:=ap;#+wx/6; Ay[j,k+3]:=bp;#+wy/6; Ax[j,k+4]:=ap-wx/6-zx/6; Ay[j,k+4]:=bp-wy/6-zy/6; Ax[j,k+5]:=ap; Ay[j,k+5]:=bp; b:=b+dy od; a:=a+dx od end: macro( INITCOND = {[constant,constant],[constant,constant,constant]} ); phaseplot:=proc(F,H,V) local i,j,p,q,r,s,phasevects,initlist,precision, pv,init,A,flow,f1,f2,npts,nsteps,m,n,step, h,cargs,carg,Alist,Ax,Ay,CC,flows: options `Copyright 1991 by Daniel Schwalbe`; if not type(F,[{procedure,numeric},{procedure,numeric}]) then ERROR(`1st argument (functions) should be a list of two procedures`) fi; if not type([H,V],[constant..constant,constant..constant]) then ERROR(`2nd and 3rd arguments (domains) must be real ranges`) fi; step:=proc(x) if x<0 then 0 else 1 fi end: nsteps:=50:npts:=1;h:=0.1:m:=16:n:=14: precision := max(Digits,trunc(evalhf(Digits))); if type(F[1],numeric) then f1:=subs(_body=F[1],proc() _body end) else f1:= subs(Heaviside=step,eval(F[1])) fi; if type(F[2],numeric) then f2:=subs(_body=F[2],proc() _body end) else f2:= subs(Heaviside=step,eval(F[2])) fi; p:=evalf(lhs(H),precision); q:=evalf(rhs(H),precision); r:=evalf(lhs(V),precision); s:=evalf(rhs(V),precision); cargs:={args[4..nargs]}; initlist:={}; for carg in cargs do if type(carg,'set') then if type(carg,'set(INITCOND)') then initlist := carg else ERROR(`initial conditions must be a set of points of type`,INITCOND) fi; elif type(carg,'list') then # A single initial condition if type(carg,INITCOND) then initlist:={carg} else ERROR(`initial condition must be of type`,INITCOND) fi; elif not type(carg,equation) then ERROR(`optional arguments must be equations`) elif op(1,carg)='stepsize' then h:=evalf(op(2,carg),precision); if not type(h,numeric) then ERROR(`stepsize must be a numerical value`) fi; elif op(1,carg)='numsteps' then nsteps:=op(2,carg); if not type(nsteps,posint) then ERROR(`numsteps must be a positive integer`) fi; elif op(1,carg)='iterations' then npts:=op(2,carg); if not type(npts,posint) then ERROR(`number of iterations must be a positive integer`) fi; elif op(1,carg)='grid' then if not type(rhs(carg),[integer,integer]) then ERROR(`grid must be a list of two positive integers`) fi; m:=op(2,carg)[1]; n:=op(2,carg)[2] else ERROR(`unknown option`,lhs(carg)) fi od; if m>0 and n>0 then Ax:=array(0..m,0..6*n+5); Ay:=array(0..m,0..6*n+5); pv:=traperror(evalhf(computevectshf(var(Ax),var(Ay), f1,f2,p,q,r,s,m,n))); if lasterror=pv then #print(`computevectshf failed`); computevects(Ax,Ay,f1,f2,p,q,r,s,m,n) fi; phasevects:=seq(seq(CURVE(BLACK,LINE, [Ax[j,6*i],Ay[j,6*i],Ax[j,6*i+1],Ay[j,6*i+1], Ax[j,6*i+2],Ay[j,6*i+2],Ax[j,6*i+3],Ay[j,6*i+3], Ax[j,6*i+4],Ay[j,6*i+4],Ax[j,6*i+5],Ay[j,6*i+5]]), j=0..m),i=0..n): else phasevects:=[0,0] fi; i:=1; for init in initlist do A:=array(-2*nsteps..2*nsteps+1); CC:=array(1..10): CC[1]:=h: if nops(init)=3 then CC[2]:=evalf(init[1]); CC[3]:=evalf(init[2]); CC[4]:=evalf(init[3]); else CC[2]:=0; CC[3]:=evalf(init[1]); CC[4]:=evalf(init[2]) fi; CC[5]:=p: CC[6]:=q: CC[7]:=r: CC[8]:=s: pv:=traperror(evalhf(phaserk( var(A),f1,f2,var(CC),nsteps,npts))); if pv = lasterror then phaserk(A,f1,f2,CC,nsteps,npts) fi; Alist:=seq(A[j],j=round(CC[9])..round(CC[10])): flow[i]:= CURVE(RED,LINE,[Alist]); i:=i+1; od; flows:=seq(flow[j],j=1..i-1): m:=max(m,16);n:=max(n,12); PLOT(AXIS(HORIZONTAL,SCALE(4,LINEAR,NUMBER)), AXIS(VERTICAL,SCALE(4,LINEAR,NUMBER)), RANGE(HORIZONTAL,p-(q-p)/m..q+(q-p)/m), RANGE(VERTICAL,r-(s-r)/n..s+(s-r)/n), phasevects,flows): end: ################################################## picard:=proc(f,init,n) local parray,ptemp,t0,y0: options `Copyright 1991 by Daniel Schwalbe`; t0:=init[1]; y0:=init[2]; ptemp:=y0+int(f('s',y0),'s'=t0..'t'); parray:=ptemp; from 2 to n do ptemp:=y0+int(f('s',subs('t'='s',ptemp)),'s'=t0..'t'); parray:=parray,sort(ptemp) od; parray; end: ################################################## #save `ODE.m`; #quit