/* Two famous ODEs and their phase portraits: Rick Wicklin, The DO Loop, Permalink: http://blogs.sas.com/content/iml/2013/09/11/create-phase-portraits-in-sas/ ?*/ /* 1. Simple harmonic oscillator: F = -kx ==> x'' = -cx -OR- Let v = x'; Then */ /* in matrix language: return( {0 1, -1 0}*z ); */ proc iml; start SHO(t, z); dxdt = z[2]; dvdt = -z[1]; return( dxdt // dvdt ); finish; x0 = {2, 0}; /* initial values */ t = do(0, 10, 0.2); /* compute solution at these time pts */ h = {1.e-6 1 1e-3}; /* min, max, and initial step size */ call ode(soln, "SHO", x0, t, h); traj = t` || (x0` // soln`); create SHO from traj[c={"Time" "x" "v"}]; append from traj; close SHO; submit; title "Simple Harmonic Motion"; proc sgplot data=SHO; label x="Position" v="Velocity"; series x=Time y=x; series x=Time y=v / y2axis; run; endsubmit; /* Phase portrait To facilitate using the GROUP= statement in PROC SGPLOT, store trajectories in "long" form as a data set with variables IC, T, X, and V. */ /* Choose initial conditions */ x0 = T(do(0.2, 3.0, 0.4)); /* initial positions */ numIC = nrow(x0); IC = x0 || j(numIC,1,0); /* each row specifies (x0, v0) */ /* Choose maximum time length of integration for each IC. You might want some trajectories to be longer than others. */ maxT = j(numIC, 1, 6.4); /* all times the same for SHO */ names = {"IC" "Time" "x" "v"}; M = {. . . .}; /* establish M as numeric matrix with 4 cols */ create PhasePortrait from M[c=names]; /* open for writing */ do i = 1 to numIC; z0 = IC[i,]`; t = do(0, maxT[i], 0.2); /* time intervals for this traj */ call ode(soln, "SHO", z0, t, h); /* fill M with data */ M = j(ncol(t),1,IC[i,1]) || t` || (z0` // soln`); append from M; /* write this trajectory */ end; close PhasePortrait; submit; /* the ASPECT= option is a SAS 9.4 feature */ title "Phase Portrait of Simple Harmonic Oscillator"; proc sgplot data=PhasePortrait aspect=1 noautolegend; series x=x y=v / group=IC; xaxis grid label="Position"; yaxis grid label="Velocity"; run; endsubmit; /***************************************/ /* 2. The Pendulum */ start Pendulum(t, z); dThetadt = z[2]; dvdt = -sin(z[1]); return( dThetadt // dvdt ); finish; x0 = T(do(0.2, 3.0, 0.4)); /* initial positions */ numIC = nrow(x0); IC = x0 || j(numIC,1,0); /* each row specifies (x0, v0) */ maxT = {6.4, 6.6, 6.8, 7.2, 8, 9, 11, 16.2}; create PhasePortrait from M[c=names]; /* open for writing */ do i = 1 to numIC; z0 = IC[i,]`; t = do(0, maxT[i], 0.2); /* time intervals for this traj */ call ode(soln, "Pendulum", z0, t, h); /* fill M with data */ M = j(ncol(t),1,IC[i,1]) || t` || (z0` // soln`); append from M; /* write this trajectory */ end; close PhasePortrait; submit; /* the ASPECT= option is a SAS 9.4 feature */ title "Phase Portrait of Pendulum"; proc sgplot data=PhasePortrait aspect=0.66 noautolegend; series x=x y=v / group=IC; xaxis grid label="Position"; yaxis grid label="Velocity"; run; quit; endsubmit;