Rzut ukośny w jednorodnym polu grawitacyjnym z uwzględnieniem siły oporu powietrza proporcjonalnej do pewnej potęgi prędkości

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

v[t_]=r'[t]

a[t_]=r''[t]

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

(* Siła grawitacji *)

Fg[t_]=m*g

(* siła oporu powietrza *)

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

(* całkowita siła *)

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

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

Przypadek γ=0 (siła oporu jest wprost proporcjonalna do prędkości), dla którego mamy analityczne rozwiązanie

γ=0

co=FullSimplify[DSolve[rownania,{x[t],y[t]},t],m>0&&g0>0&&k>0]

(* Dostaliśmy tzw. całkę ogólną układu równań różniczkowych, zależną od 4 dowolnych stałych, C[1], C[2], ..., C[4]. Te stałe wyliczymy, zakładając tzw. warunki początkowe: wektor położenia w chwili t=0 i wektor prędkości w chwili t=0.
Te warunki początkowe formułujemy zupełnie ogólnie:
Położenie w chwili początkowej: r[0]={x0,y0}.
Prędkość w chwili początkowej: v[0]={vx0,vy0}.
 *)

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

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

(* To już jest całka szczególna układu równańn różniczkowych, gdzie stałe C[1], C[2], ..., C[4] zostały już wyeliminowane.
*)

(* Następne komendy pokazuja, jak z rozwiazania (cs) uzyskac potrzebne nam wzory *)

cs[[1]][[1]]

cs[[1]][[2]]

x1[t_]=x[t]/.cs[[1]][[1]]

y1[t_]=y[t]/.cs[[1]][[2]]

(* Przykladowe dane dla rysunku *)

θ=Pi/3;

m=0.5; (* kg *)

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

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

v0=20; (* m/s *)

h=50; (* m *)

(* dla porownania znane rozwiazanie bez oporu powietrza *)

xsw[t_]=t*v0* Cos[θ]

ysw[t_]=h-(1/2)*g0*t*t+t*v0* Sin[θ]

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

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

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

(* Przykładowe dane dla rysunku *)

γ=1.2;

θ=Pi/3;

m=0.5; (* kg *)

k=0.01; (* N*(s/m)^(γ+1) *)

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

v0=20; (* m/s *)

h=50; (* m *)

ClearAll[x,y]

rownania

WP

num=NDSolve[{rownania,WP},{x,y},{t,0,15}];

num[[1]][[1]]

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

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

ParametricPlot[{{xsw[t],ysw[t]},{x1[t],y1[t]},{x2[t],y2[t]}},{t,0,6},PlotRange->{{0,60},{-5,70}},AxesLabel->{"x","y"},PlotLegends->"Expressions"]