Example 2.1: Free fall (Ch 2, page 25)

Clear["Global`*"]
Off[General::spell1]

<<VectorFieldPlots`;

System

Parameters

{(Subscript[x, min]),(Subscript[x, max])}={-4,4};
{(Subscript[y, min]),(Subscript[y, max])}={-4,4};

f[x_,y_]:=y

g[x_,y_]:=-9.81

Vector field

V={f[x,y],g[x,y]}

Jacobian

a11=Simplify[D[f[x,y],x]]

a12=Simplify[D[f[x,y],y]]

a21=Simplify[D[g[x,y],x]]

a22=Simplify[D[g[x,y],y]]

J={{a11,a12},{a21,a22}}

MatrixForm[J]

Acceleration

γ=J . V

Phase Portrait

x=.
y=.
t=.

(* This procedure enables to reset the values attributed at x, y and t intableplot *)

k = 1 ; (x(0), y(0)) = (0.8, 4)

solution=NDSolve[{x'[t]==f[x[t],y[t]],x[0]==0.8,

y'[t]==g[x[t],y[t]],y[0]==4},{x,y},

{t,0,100},Method->"ExplicitRungeKutta", MaxSteps->Infinity,StartingStepSize->0.001];phaseplot0=ParametricPlot[Evaluate[{x[t],y[t]}/.solution],{t,0,100},PlotPoints->1000,PlotRange->{{(Subscript[x, min]),(Subscript[x, max])},{(Subscript[y, min]),(Subscript[y, max])}},AxesLabel->{"X","Y"},PlotStyle->{{Red,Thickness[0.005]}},AspectRatio->1,Frame->True];

tableplot={};
For[t=0,t≤5,t=t+ .05,
x=Evaluate[{x[t],y[t]}/.solution][[1]][[1]];
y=Evaluate[{x[t],y[t]}/.solution][[1]][[2]];
vx=f[x,y];
vy=g[x,y];
tmptable={{x,y},{vx,vy}};
tableplot=Append[tableplot,tmptable];];
flow0=ListVectorFieldPlot[tableplot ,ScaleFunction->(.25&),ColorFunction->(Hue[0.66]&),Frame->True,AspectRatio->1]

k = 2 ; (x(0), y(0)) = (1.6, 4)

x=.
y=.
t=.

solution=NDSolve[{x'[t]==f[x[t],y[t]],x[0]==1.6,

y'[t]==g[x[t],y[t]],y[0]==4},{x,y},

{t,0,100},Method->"ExplicitRungeKutta", MaxSteps->Infinity,StartingStepSize->0.001];phaseplot1=ParametricPlot[Evaluate[{x[t],y[t]}/.solution],{t,0,100},PlotPoints->1000,PlotRange->{{(Subscript[x, min]),(Subscript[x, max])},{(Subscript[y, min]),(Subscript[y, max])}},AxesLabel->{"X","Y"},PlotStyle->{{Red,Thickness[0.005]}},AspectRatio->1,Frame->True];

tableplot={};
For[t=0,t≤5,t=t+ .05,
x=Evaluate[{x[t],y[t]}/.solution][[1]][[1]];
y=Evaluate[{x[t],y[t]}/.solution][[1]][[2]];
vx=f[x,y];
vy=g[x,y];
tmptable={{x,y},{vx,vy}};
tableplot=Append[tableplot,tmptable];];
flow1=ListVectorFieldPlot[tableplot ,ScaleFunction->(.25&),ColorFunction->(Hue[0.66]&),Frame->True,AspectRatio->1]

k = 0 ; (x(0), y(0)) = (0, 4)

x=.
y=.
t=.

solution=NDSolve[{x'[t]==f[x[t],y[t]],x[0]==0,

y'[t]==g[x[t],y[t]],y[0]==4},{x,y},

{t,0,100},Method->"ExplicitRungeKutta", MaxSteps->Infinity,StartingStepSize->0.001];phaseplot2=ParametricPlot[Evaluate[{x[t],y[t]}/.solution],{t,0,100},PlotPoints->1000,PlotRange->{{(Subscript[x, min]),(Subscript[x, max])},{(Subscript[y, min]),(Subscript[y, max])}},AxesLabel->{"X","Y"},PlotStyle->{{Red,Thickness[0.005]}},AspectRatio->1,Frame->True];

tableplot={};
For[t=0,t≤5,t=t+ .05,
x=Evaluate[{x[t],y[t]}/.solution][[1]][[1]];
y=Evaluate[{x[t],y[t]}/.solution][[1]][[2]];
vx=f[x,y];
vy=g[x,y];
tmptable={{x,y},{vx,vy}};
tableplot=Append[tableplot,tmptable];];
flow2=ListVectorFieldPlot[tableplot ,ScaleFunction->(.25&),ColorFunction->(Hue[0.66]&),Frame->True,AspectRatio->1]

k = – 1 ; (x(0), y(0)) = (– 0.8, 4)

x=.
y=.
t=.

solution=NDSolve[{x'[t]==f[x[t],y[t]],x[0]==-0.8,

y'[t]==g[x[t],y[t]],y[0]==4},{x,y},

{t,0,100},Method->"ExplicitRungeKutta", MaxSteps->Infinity,StartingStepSize->0.001];phaseplot3=ParametricPlot[Evaluate[{x[t],y[t]}/.solution],{t,0,100},PlotPoints->1000,PlotRange->{{(Subscript[x, min]),(Subscript[x, max])},{(Subscript[y, min]),(Subscript[y, max])}},AxesLabel->{"X","Y"},PlotStyle->{{Red,Thickness[0.005]}},AspectRatio->1,Frame->True];

tableplot={};
For[t=0,t≤5,t=t+ .05,
x=Evaluate[{x[t],y[t]}/.solution][[1]][[1]];
y=Evaluate[{x[t],y[t]}/.solution][[1]][[2]];
vx=f[x,y];
vy=g[x,y];
tmptable={{x,y},{vx,vy}};
tableplot=Append[tableplot,tmptable];];
flow3=ListVectorFieldPlot[tableplot ,ScaleFunction->(.25&),ColorFunction->(Hue[0.66]&),Frame->True,AspectRatio->1]

k = – 2 ; (x(0), y(0)) = (– 1.6, 4)

x=.
y=.
t=.

solution=NDSolve[{x'[t]==f[x[t],y[t]],x[0]==-1.6,

y'[t]==g[x[t],y[t]],y[0]==4},{x,y},

{t,0,100},Method->"ExplicitRungeKutta", MaxSteps->Infinity,StartingStepSize->0.001];phaseplot4=ParametricPlot[Evaluate[{x[t],y[t]}/.solution],{t,0,100},PlotPoints->1000,PlotRange->{{(Subscript[x, min]),(Subscript[x, max])},{(Subscript[y, min]),(Subscript[y, max])}},AxesLabel->{"X","Y"},PlotStyle->{{Red,Thickness[0.005]}},AspectRatio->1,Frame->True];

tableplot={};
For[t=0,t≤5,t=t+ .05,
x=Evaluate[{x[t],y[t]}/.solution][[1]][[1]];
y=Evaluate[{x[t],y[t]}/.solution][[1]][[2]];
vx=f[x,y];
vy=g[x,y];
tmptable={{x,y},{vx,vy}};
tableplot=Append[tableplot,tmptable];];
flow4=ListVectorFieldPlot[tableplot ,ScaleFunction->(.25&),ColorFunction->(Hue[0.66]&),Frame->True,AspectRatio->1]

k = – 3 ; (x(0), y(0)) = (– 2.4, 4)

x=.
y=.
t=.

solution=NDSolve[{x'[t]==f[x[t],y[t]],x[0]==-2.4,

y'[t]==g[x[t],y[t]],y[0]==4},{x,y},

{t,0,100},Method->"ExplicitRungeKutta", MaxSteps->Infinity,StartingStepSize->0.001];phaseplot5=ParametricPlot[Evaluate[{x[t],y[t]}/.solution],{t,0,100},PlotPoints->1000,PlotRange->{{(Subscript[x, min]),(Subscript[x, max])},{(Subscript[y, min]),(Subscript[y, max])}},AxesLabel->{"X","Y"},PlotStyle->{{Red,Thickness[0.005]}},AspectRatio->1,Frame->True];

tableplot={};
For[t=0,t≤5,t=t+ .05,
x=Evaluate[{x[t],y[t]}/.solution][[1]][[1]];
y=Evaluate[{x[t],y[t]}/.solution][[1]][[2]];
vx=f[x,y];
vy=g[x,y];
tmptable={{x,y},{vx,vy}};
tableplot=Append[tableplot,tmptable];];
flow5=ListVectorFieldPlot[tableplot ,ScaleFunction->(.25&),ColorFunction->(Hue[0.66]&),Frame->True,AspectRatio->1]

k = – 4 ; (x(0), y(0)) = (– 4.2, 4)

x=.
y=.
t=.

solution=NDSolve[{x'[t]==f[x[t],y[t]],x[0]==-3.2,

y'[t]==g[x[t],y[t]],y[0]==4},{x,y},

{t,0,100},Method->"ExplicitRungeKutta", MaxSteps->Infinity,StartingStepSize->0.001];phaseplot6=ParametricPlot[Evaluate[{x[t],y[t]}/.solution],{t,0,100},PlotPoints->1000,PlotRange->{{(Subscript[x, min]),(Subscript[x, max])},{(Subscript[y, min]),(Subscript[y, max])}},AxesLabel->{"X","Y"},PlotStyle->{{Red,Thickness[0.005]}},AspectRatio->1,Frame->True];

tableplot={};
For[t=0,t≤5,t=t+ .05,
x=Evaluate[{x[t],y[t]}/.solution][[1]][[1]];
y=Evaluate[{x[t],y[t]}/.solution][[1]][[2]];
vx=f[x,y];
vy=g[x,y];
tmptable={{x,y},{vx,vy}};
tableplot=Append[tableplot,tmptable];];
flow6=ListVectorFieldPlot[tableplot ,ScaleFunction->(.25&),ColorFunction->(Hue[0.66]&),Frame->True,AspectRatio->1]

Show[phaseplot0,flow0,phaseplot1,flow1,phaseplot2,flow2,phaseplot3,flow3,phaseplot4,flow4,phaseplot5,flow5,phaseplot6,flow6,ImageSize->{600,600}]