Glider

Case 1: D=0 (No air drag)

PhasePortrait

eps=.01;
StreamPlot[{1/v*(-Cos[θ]+v^2),-Sin[θ]},{θ,-π,π},{v,eps,3},FrameLabel->{"θ","v"}]

Solution Curves on Cylinder

plota=ParametricPlot3D[{3/2*Cos[θ],3/2*Sin[θ],v},{θ,-π,π},{v,eps,3},Boxed->False,Axes->False,Mesh->False,PlotStyle->Opacity[0.6],PlotPoints->20];
For[i=eps,i≤6+eps,i++,
v0=i/2;
s=NDSolve[{v'[t]==-Sin[θ[t]],θ'[t]==1/v[t]*(-Cos[θ[t]]+v[t]^2),θ[0]==0,v[0]==v0},{v,θ},{t,0,2*π}];
plotb=ParametricPlot3D[Evaluate[{(3/2+.001)*Cos[θ[t]],(3/2+.001)*Sin[θ[t]],v[t]}/.s],{t,0,2*π},PlotStyle->{Thick,Red}];
plota=Show[plota,plotb]];
plotc=ParametricPlot3D[{(3/2+.001)*Cos[θ],(3/2+.001)*Sin[θ],0},{θ,0,2*π},PlotStyle->{Thick,Blue}];
plotd=ParametricPlot3D[{(3/2+.001)*Cos[θ],(3/2+.001)*Sin[θ],Sqrt[Cos[θ]]},{θ,0,2*π},PlotStyle->{Thick,Green}];
plotf=ParametricPlot3D[{(3/2+.001)*Cos[0],(3/2+.001)*Sin[0],v},{v,eps,3},PlotStyle->{Thick,Purple}];
plotg=ParametricPlot3D[{(3/2+.001)*Cos[π],(3/2+.001)*Sin[π],v},{v,eps,3},PlotStyle->{Thick,Purple}];
Show[plota,plotc,plotd,plotf,plotg,ImageSize->Large,ViewPoint->{6,3,2}]

Plot of Glider Dynamics

Since ((((dx/dt)))^2)+((((dy/dt)))^2)=((((dv/dt)))^2)we can recover the dynamics of the actual glider. Namely, we obtain the 4-d system 
(dx/dt)=v(t)cos(θ(t))
(dy/dt)=v(t)sin(θ(t))
(dv/dt)=-sin(θ(t))-(Dv^2)
v((dθ)/dt)=-Cos(θ(t))+(v^2)

We can use the phase portrait to give us interesting initial conditions for this system.

v0=2;
θ0=0;
s1=NDSolve[{x'[t]==v[t]*Cos[θ[t]],y'[t]==v[t]*Sin[θ[t]],v'[t]==-Sin[θ[t]],θ'[t]==1/v[t]*(-Cos[θ[t]]+v[t]^2),x[0]==0,y[0]==0,θ[0]==θ0,v[0]==v0},{x,y,v,θ},{t,0,20}];
plot1=ParametricPlot[Evaluate[{x[t],y[t]}/.s1],{t,0,20},PlotRange->{{0,20},{-2,2}},PlotStyle->{Thick,Red}];
v0=1;
θ0=0;
s1=NDSolve[{x'[t]==v[t]*Cos[θ[t]],y'[t]==v[t]*Sin[θ[t]],v'[t]==-Sin[θ[t]],θ'[t]==1/v[t]*(-Cos[θ[t]]+v[t]^2),x[0]==0,y[0]==0,θ[0]==θ0,v[0]==v0},{x,y,v,θ},{t,0,20}];
plot2=ParametricPlot[Evaluate[{x[t],y[t]}/.s1],{t,0,20},PlotRange->{{0,20},{-2,2}},PlotStyle->{Thick,Green}];
v0=1/2;
θ0=0;
s1=NDSolve[{x'[t]==v[t]*Cos[θ[t]],y'[t]==v[t]*Sin[θ[t]],v'[t]==-Sin[θ[t]],θ'[t]==1/v[t]*(-Cos[θ[t]]+v[t]^2),x[0]==0,y[0]==0,θ[0]==θ0,v[0]==v0},{x,y,v,θ},{t,0,20}];
plot3=ParametricPlot[Evaluate[{x[t],y[t]}/.s1],{t,0,20},PlotRange->{{0,20},{-2,2}},PlotStyle->{Thick,Blue}];
Show[plot1,plot2,plot3]

It is clear that the dynamics of the glider is heavily influenced by the intitial velocity. It the velocity is greater than 1 we waste a lot of time flying backwards while if the velocity is less than one we we spend a lot of time “dipping”.

Lets now explore some extreme initial conditions.

v0=4;
θ0=0;
s1=NDSolve[{x'[t]==v[t]*Cos[θ[t]],y'[t]==v[t]*Sin[θ[t]],v'[t]==-Sin[θ[t]],θ'[t]==1/v[t]*(-Cos[θ[t]]+v[t]^2),x[0]==0,y[0]==0,θ[0]==θ0,v[0]==v0},{x,y,v,θ},{t,0,20}];
plot1=ParametricPlot[Evaluate[{x[t],y[t]}/.s1],{t,0,20},PlotRange->{{0,20},{-2,2}},PlotStyle->{Thick,Red}];
v0=1;
θ0=0;
s1=NDSolve[{x'[t]==v[t]*Cos[θ[t]],y'[t]==v[t]*Sin[θ[t]],v'[t]==-Sin[θ[t]],θ'[t]==1/v[t]*(-Cos[θ[t]]+v[t]^2),x[0]==0,y[0]==0,θ[0]==θ0,v[0]==v0},{x,y,v,θ},{t,0,20}];
plot2=ParametricPlot[Evaluate[{x[t],y[t]}/.s1],{t,0,20},PlotRange->{{0,20},{-2,2}},PlotStyle->{Thick,Green}];
v0=1/100;
θ0=0;
s1=NDSolve[{x'[t]==v[t]*Cos[θ[t]],y'[t]==v[t]*Sin[θ[t]],v'[t]==-Sin[θ[t]],θ'[t]==1/v[t]*(-Cos[θ[t]]+v[t]^2),x[0]==0,y[0]==0,θ[0]==θ0,v[0]==v0},{x,y,v,θ},{t,0,20}];
plot3=ParametricPlot[Evaluate[{x[t],y[t]}/.s1],{t,0,20},PlotRange->{{0,20},{-2,2}},PlotStyle->{Thick,Blue}];
Show[plot1,plot2,plot3]

This is awesome!!! It looks like the trajectories become circular if the velocity is too large. If the velocity is too small it seems almost like the trajectories are cycloids. Can this be proved?? The circular case makes sense if (Subscript[v, 0])>>1 then (dv/(dθ))≈0 which implies that (dx/dt)≈(Subscript[v, 0])Cos(θ(t)), (dy/dt)≈(Subscript[v, 0])Sin(θ(t)), i.e. circular motion. What about if (Subscript[v, 0])<<1. This would require analysis of the heteroclinic orbit connecting the two points on the curve v=0 that act like saddle points.

Case 2: D>0 (Adding Air Drag)

PhasePortrait

d=.1;
eps=.01;
StreamPlot[{1/v*(-Cos[θ]+v^2),-Sin[θ]-d*v^2},{θ,-π,π},{v,eps,3},FrameLabel->{"θ","v"}]

Solution Curves on Cylinder

plota=ParametricPlot3D[{3/2*Cos[θ],3/2*Sin[θ],v},{θ,-π,π},{v,eps,3},Boxed->False,Axes->False,Mesh->False,PlotStyle->Opacity[0.6],PlotPoints->20];
For[i=eps,i≤6+eps,i++,
v0=i/2;
s=NDSolve[{v'[t]==-Sin[θ[t]]-d*v[t]^2,θ'[t]==1/v[t]*(-Cos[θ[t]]+v[t]^2),θ[0]==0,v[0]==v0},{v,θ},{t,0,20}];
plotb=ParametricPlot3D[Evaluate[{(3/2+.001)*Cos[θ[t]],(3/2+.001)*Sin[θ[t]],v[t]}/.s],{t,0,20},PlotStyle->{Thick,Red}];
plota=Show[plota,plotb]];
plotc=ParametricPlot3D[{(3/2+.001)*Cos[θ],(3/2+.001)*Sin[θ],0},{θ,0,2*π},PlotStyle->{Thick,Blue}];
plotd=ParametricPlot3D[{(3/2+.001)*Cos[θ],(3/2+.001)*Sin[θ],Sqrt[Cos[θ]]},{θ,0,2*π},PlotStyle->{Thick,Green}];
plotf=ParametricPlot3D[{(3/2+.001)*Cos[θ],(3/2+.001)*Sin[θ],Sqrt[-Sin[θ]/d]},{θ,-π,π},PlotStyle->{Thick,Purple}];
Show[plota,plotc,plotd,plotf,ImageSize->Large,ViewPoint->{6,3,2}]

Plot of Glider Dynamics

Since ((((dx/dt)))^2)+((((dy/dt)))^2)=((((dv/dt)))^2)we can recover the dynamics of the actual glider. Namely, we obtain the 4-d system 
(dx/dt)=v(t)cos(θ(t))
(dy/dt)=v(t)sin(θ(t))
(dv/dt)=-sin(θ(t))-(Dv^2)
v((dθ)/dt)=-Cos(θ(t))+(v^2)

We can use the phase portrait to give us interesting initial conditions for this system.

d=.1;
StreamPlot[{1/v*(-Cos[θ]+v^2),-Sin[θ]-d*v^2},{θ,-π,π},{v,-3,3},FrameLabel->{"θ","v"}]
v0=2;
θ0=0;
s1=NDSolve[{x'[t]==v[t]*Cos[θ[t]],y'[t]==v[t]*Sin[θ[t]],v'[t]==-Sin[θ[t]]-d*v[t]^2,θ'[t]==1/v[t]*(-Cos[θ[t]]+v[t]^2),x[0]==0,y[0]==0,θ[0]==θ0,v[0]==v0},{x,y,v,θ},{t,0,20}];
plot1=ParametricPlot[Evaluate[{x[t],y[t]}/.s1],{t,0,20},PlotRange->{{0,20},{-2,2}},PlotStyle->{Thick,Red}];
v0=1;
θ0=0;
s1=NDSolve[{x'[t]==v[t]*Cos[θ[t]],y'[t]==v[t]*Sin[θ[t]],v'[t]==-Sin[θ[t]]-d*v[t]^2,θ'[t]==1/v[t]*(-Cos[θ[t]]+v[t]^2),x[0]==0,y[0]==0,θ[0]==θ0,v[0]==v0},{x,y,v,θ},{t,0,20}];
plot2=ParametricPlot[Evaluate[{x[t],y[t]}/.s1],{t,0,20},PlotRange->{{0,20},{-2,2}},PlotStyle->{Thick,Green}];
v0=1/2;
θ0=0;
s1=NDSolve[{x'[t]==v[t]*Cos[θ[t]],y'[t]==v[t]*Sin[θ[t]],v'[t]==-Sin[θ[t]]-d*v[t]^2,θ'[t]==1/v[t]*(-Cos[θ[t]]+v[t]^2),x[0]==0,y[0]==0,θ[0]==θ0,v[0]==v0},{x,y,v,θ},{t,0,20}];
plot3=ParametricPlot[Evaluate[{x[t],y[t]}/.s1],{t,0,20},PlotRange->{{0,20},{-2,2}},PlotStyle->{Thick,Blue}];
Show[plot1,plot2,plot3]

Now it is clear that a little bit of drag causes the glider to begin descending.

d=.01;
StreamPlot[{1/v*(-Cos[θ]+v^2),-Sin[θ]-d*v^2},{θ,-π,π},{v,-3,3},FrameLabel->{"θ","v"}]
v0=2;
θ0=0;
s1=NDSolve[{x'[t]==v[t]*Cos[θ[t]],y'[t]==v[t]*Sin[θ[t]],v'[t]==-Sin[θ[t]]-d*v[t]^2,θ'[t]==1/v[t]*(-Cos[θ[t]]+v[t]^2),x[0]==0,y[0]==0,θ[0]==θ0,v[0]==v0},{x,y,v,θ},{t,0,20}];
plot1=ParametricPlot[Evaluate[{x[t],y[t]}/.s1],{t,0,20},PlotRange->{{0,20},{-2,2}},PlotStyle->{Thick,Red}];
v0=1;
θ0=0;
s1=NDSolve[{x'[t]==v[t]*Cos[θ[t]],y'[t]==v[t]*Sin[θ[t]],v'[t]==-Sin[θ[t]]-d*v[t]^2,θ'[t]==1/v[t]*(-Cos[θ[t]]+v[t]^2),x[0]==0,y[0]==0,θ[0]==θ0,v[0]==v0},{x,y,v,θ},{t,0,20}];
plot2=ParametricPlot[Evaluate[{x[t],y[t]}/.s1],{t,0,20},PlotRange->{{0,20},{-2,2}},PlotStyle->{Thick,Green}];
v0=1/2;
θ0=0;
s1=NDSolve[{x'[t]==v[t]*Cos[θ[t]],y'[t]==v[t]*Sin[θ[t]],v'[t]==-Sin[θ[t]]-d*v[t]^2,θ'[t]==1/v[t]*(-Cos[θ[t]]+v[t]^2),x[0]==0,y[0]==0,θ[0]==θ0,v[0]==v0},{x,y,v,θ},{t,0,20}];
plot3=ParametricPlot[Evaluate[{x[t],y[t]}/.s1],{t,0,20},PlotRange->{{0,20},{-2,2}},PlotStyle->{Thick,Blue}];
Show[plot1,plot2,plot3]

d=.5;
StreamPlot[{1/v*(-Cos[θ]+v^2),-Sin[θ]-d*v^2},{θ,-π,π},{v,-3,3},FrameLabel->{"θ","v"}]
v0=2;
θ0=0;
s1=NDSolve[{x'[t]==v[t]*Cos[θ[t]],y'[t]==v[t]*Sin[θ[t]],v'[t]==-Sin[θ[t]]-d*v[t]^2,θ'[t]==1/v[t]*(-Cos[θ[t]]+v[t]^2),x[0]==0,y[0]==0,θ[0]==θ0,v[0]==v0},{x,y,v,θ},{t,0,20}];
plot1=ParametricPlot[Evaluate[{x[t],y[t]}/.s1],{t,0,20},PlotRange->{{0,20},{-2,2}},PlotStyle->{Thick,Red}];
v0=1;
θ0=0;
s1=NDSolve[{x'[t]==v[t]*Cos[θ[t]],y'[t]==v[t]*Sin[θ[t]],v'[t]==-Sin[θ[t]]-d*v[t]^2,θ'[t]==1/v[t]*(-Cos[θ[t]]+v[t]^2),x[0]==0,y[0]==0,θ[0]==θ0,v[0]==v0},{x,y,v,θ},{t,0,20}];
plot2=ParametricPlot[Evaluate[{x[t],y[t]}/.s1],{t,0,20},PlotRange->{{0,20},{-2,2}},PlotStyle->{Thick,Green}];
v0=1/2;
θ0=0;
s1=NDSolve[{x'[t]==v[t]*Cos[θ[t]],y'[t]==v[t]*Sin[θ[t]],v'[t]==-Sin[θ[t]]-d*v[t]^2,θ'[t]==1/v[t]*(-Cos[θ[t]]+v[t]^2),x[0]==0,y[0]==0,θ[0]==θ0,v[0]==v0},{x,y,v,θ},{t,0,20}];
plot3=ParametricPlot[Evaluate[{x[t],y[t]}/.s1],{t,0,20},PlotRange->{{0,20},{-2,2}},PlotStyle->{Thick,Blue}];
Show[plot1,plot2,plot3]

How does superman remain horizontal while he is flying????? Clearly superman must be able to eliminate gravity and drag around him. This is much more powerful than just being super strong!! If he could do that why doesn’t he could just launch people into outerspace by being near them.