Rzut ukośny w jednorodnym ziemskim polu grawitacyjnym z uwzględnieniem siły oporu powietrza proporcjonalnej do pewnej potęgi prędkości. Bez straty ogólności możemy założyć, że ruch odbywa się w płaszczyźnie xy

(* wektor położenia *)

r[t_]={x[t],y[t]}

(* prędkość *)

v[t_]=r'[t]

(* przyspieszenie *)

a[t_]=r''[t]

(* wektor przyspieszenia ziemskiego *)

g={0,-g0}  (* g0 > 0 ! *)

(* Siła grawitacji *)

Fg[t_]=m*g

(* Siła oporu powietrza jest skierowana przeciwnie do wektora prędkości, a jej długość wynosi k*|v[t]|^(γ+1)  *)

Fop[t_]=-k*(v[t].v[t])^(γ/2)*v[t]

(* całkowita siła *)

F[t_]=Fg[t]+Fop[t]

(* II zasada dynamiki Newtona *)

rownania=m*a[t]==F[t]

Interesuje nas rozwiązanie z konkretnymi warunkami początkowymi: podajemy wektor położenia w chwili t=0 i wektor prędkości w chwili t=0.

WP={x[0]==0,y[0]==h,x'[0]==v0*Cos[θ],y'[0]==v0*Sin[θ]}

Przypadek bez oporu powietrza (k=0), dla którego mamy najprostsze analityczne rozwiązanie

cs0=FullSimplify[DSolve[{rownania,WP}/.{k->0},{x[t],y[t]},t],m>0&&g0>0&&k≥ 0&&v0>0&&Pi/2>θ>0]

(* Następne komendy pokazują, jak z rozwiązania (cs0) uzyskać potrzebne nam wzory *)

cs0[[1]][[1]]

cs0[[1]][[2]]

xsw[t_]=x[t]/.cs0[[1]][[1]]

ysw[t_]=y[t]/.cs0[[1]][[2]]

(* Przykładowe dane dla rysunku *)

θ=Pi/3;

m=0.5; (* kg *)

g0=10; (* m/s^2 *)

v0=20; (* m/s *)

h=50; (* m *)

ParametricPlot[{{xsw[t],ysw[t]}},{t,0,5.5},PlotRange->All,AxesLabel->{"x","y"}]

Przypadek z siłą oporu powietrza proporcjonalną do wektora prędkości (γ=0). Także dla tego przypadku mamy analityczne rozwiązanie.

ClearAll[x,y]

ClearAll[θ,m,g0,v0,h];

k1=0.01; (* N*(s/m) *)

cs1=FullSimplify[DSolve[{rownania,WP}/.{γ-> 0},{x[t],y[t]},t],m>0&&g0>0&&k≥ 0&&v0>0&&Pi/2>θ>0]

(* Z rozwiązania (cs1) wydobywam wzory *)

x1[t_]=(x[t]/.cs1[[1]][[1]])/.{k->k1}

y1[t_]=(y[t]/.cs1[[1]][[2]])/.{k->k1}

(* Przykladowe dane dla rysunku *)

θ=Pi/3;

m=0.5; (* kg *)

g0=10; (* m/s^2 *)

v0=20; (* m/s *)

h=50; (* m *)

ParametricPlot[{{xsw[t],ysw[t]},{x1[t],y1[t]}},{t,0,5.5},PlotRange->All,AxesLabel->{"x","y"}]

Przypadek γ≠0, dla którego nie mamy analitycznego rozwiązania i równanie różniczkowe musimy rozwiązywać numerycznie

ClearAll[x,y]

k2=0.002; (* N*(s/m)^(γ+1) *)

cs2num=NDSolve[{rownania,WP}/.{k->k2,γ->6/5},{x,y},{t,0,15}]

cs2num[[1]][[1]]

x2[t_]=Evaluate[x[t]/.cs2num[[1]][[1]]]

y2[t_]=Evaluate[y[t]/.cs2num[[1]][[2]]]

ParametricPlot[{{xsw[t],ysw[t]},{x1[t],y1[t]},{x2[t],y2[t]}},{t,0,5.5},PlotRange->All,AxesLabel->{"x","y"}]