# Kapitel 3_c: 
# Der schiefe Wurf mit Widerstand (Viskose Reibung)
# © 2015 Friedrich U. Mathiak,
# mathiak@mechanik-info.de
# 
> restart: with(LinearAlgebra):
# Beispiel 4-6:
# 
# 
> ode1:=diff(v(alpha),alpha)*cos(alpha)-v(alpha)*(sin(alpha)+k1*v(alpha)/(m*g));
> AB1:=v(alpha0)=v0;
> loe1:=dsolve({ode1,AB1});assign(%)
;
> v:=unapply(simplify(v(alpha)),alpha,k1);
# Wir überprüfen den Anfangswert:
> v(alpha0,k1);
# Für ±  gegen -À/2  geht v(±) gegen  und damit gegen die Endgeschwindigkeit
> vinf:=unapply(limit(v(alpha,k1),alpha=-Pi/2),k1);
# Im Scheitelpunkt der Wurfbahn (± = 0) hat der Geschwindigkeitsvektor nur die Horizontalkomponente
> vH:=v(0,k1);
# Fehlt mit k1 = 0 der Reibungswiderstand, dann ist mit
> v(alpha,0)*cos(alpha);
# die Horizontalkomponente der Geschwindigkeit konstant. Die Berechnung der Bahnkurve erfolgt mit x10 = x20 = 0:
> x1:=-1/g*int(v(alpha,k1)^2,alpha=alpha0..alpha);
> x2:=1/g*int(v(x,k1)^2*tan(x),x=alpha..alpha0,'AllSolutions') assuming alpha0 >= alpha;
# Für eine grafische Ausgabe der Geschwindigkeit und der Bahnkurve konkretisieren wir die Aufgabe und wählen
> m:=0.1; alpha0:=Pi/4; g:=9.81; v0:=20.; k1:=1;
> plot([v(alpha,1),v(alpha,0.5)],alpha=-Pi/2 .. Pi/4,title = "\nGeschwindigkeit eines Körpers bei viskoser Reibung\n", legend = ["k1 = 1.0", "k1 = 0.5"], titlefont = ["ARIAL", 15], labels = ["alpha", "v/v0"], labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 10], linestyle = [solid, longdash], axesfont = ["HELVETICA", "ROMAN", 8], legendstyle = [font = ["HELVETICA", 9], location = bottom],gridlines=true,axes=boxed);
> plot([x1,x2],alpha=-Pi/2 .. Pi/4,title = "\nBahnkurve eines Körpers bei viskoser Reibung (k1 = 1)\n", legend = ["x1", "x2"], titlefont = ["ARIAL", 15], labels = ["alpha", "Bahnkurve"], labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 10], linestyle = [solid, longdash], axesfont = ["HELVETICA", "ROMAN", 8], legendstyle = [font = ["HELVETICA", 9], location = bottom],gridlines=true,axes=boxed);
> m:='m': g:='g': k1:='k1': x1:='x1': x2:='x2': 
# Statt des Parameters ±, der die Neigung des Geschwindigkeitsvektors bezüglich der positiven x1-Achse beschreibt, können wir auch die Zeit t als unabhängige Variable einführen. Aus dem Schwerpunktsatz erhalten wir eine Vektorgleichung, die in zwei entkoppelte Differenzialgleichungen zerfällt. Wir berechnen zunächst die Bewegung der Masse m in x1-Richtung. Die zugehörige Differenzialgleichung ist
> ode1:=m*diff(x1(t),t,t)+k1*diff(x1(t),t);
# Die Anfangsbedingungen sind:
> AB1:=x1(0)=0,D(x1)(0)=v10;
> loe1:=dsolve({ode1,AB1},x1(t));
> x1:=unapply(rhs(loe1),t,k1);
> v1:=unapply(diff(x1(t,k1),t),t);
> a1:=unapply(diff(v1(t),t),t);
# Im nächsten Schritt wird die Bewegung in x2-Richtung berechnet. 
> ode2:=m*diff(x2(t),t,t)+k1*diff(x2(t),t) + m*g;
# Für diese Richtung sind folgende Anfangsbedingungen zu erfüllen:
> AB2:=x2(0)=0,D(x2)(0)=v20;
> loe2:=dsolve({ode2,AB2});
> x2:=unapply(rhs(loe2),t,k1);
> v2:=unapply(diff(x2(t,k1),t),t);
> a2:=unapply(diff(v2(t),t),t);
# Den Betrag des Geschwindigkeitsvektors errechnen wir zu
> v:=unapply(sqrt(v1(t)^2+v2(t)^2),t,k1);
# und für den Betrag  Beschleunigungsvektors  folgt
> a:=unapply(sqrt(a1(t)^2+a2(t)^2),t,k1);
# Die Wurfbahn ist  keine quadratische Parabel wie beim Wurf ohne Luftwiderstand. Die Scheitelkoordinaten der Wurfbahn berechnen wir wie folgt:
# Zunächst wird aus v2 = 0 die Steigzeit tH berechnet:
> tH:=solve(v2(t),t);
# Setzen wir diesen Wert in die Wurfbahnkoordinaten x1 und x2  ein, dann erhalten wir die Scheitelkoordinaten des schiefen Wurfes mit viskoser Reibung:
> x1H:=unapply(simplify(x1(tH,k1)),k1); 
> x2H:=unapply(simplify(x2(tH,k1)),k1);
# Die Wurfbahn besitzt eine vertikale Asymptote (x10 = 0), bei 
> x1max:=unapply(m/k1*v10,k1);
# und der Körper hat im Grenzfall nur die Geschwindigkeitskomponente
> v2grenz:=-m*g/k1;
# Wir führen den Parameter
> alpha:=unapply(arctan(v2(t)/v1(t)),t,k1);
# ein und konkretisieren die Aufgabe.  Die Masse m wird zum Zeitpunkt t = 0 unter dem Winkel ±0 mit der Geschwindigkeit v0 bei x10 = x20 = 0 abgeworfen. Für die nachfolgende  Zahlenrechnung wählen wir die Werte:
> m:=0.5; g:=9.81;  alpha0:=Pi/4.; v0:=20.; v10:=v0*cos(alpha0); v20:=v0*sin(alpha0); #k1:=1.;
# Die Masse m erreicht mit x2 = 0  den Boden zum Zeitpunkt 
> tW1:=solve(x2(t,0.1),t, useassumptions) assuming t>0;
> tW2:=solve(x2(t,0.5),t, useassumptions) assuming t>0;
> tW3:=solve(x2(t,1.0),t, useassumptions) assuming t>0;
# Wir geben den Verlauf der Tangentenneigung für den Widerstandsbeiwert k1 = 1.0  grafisch aus:
> plot(alpha(t,1),t=0..tW3,title = "\nNeigung der Tangente an die Bahnkurve\n",titlefont = ["ARIAL", 15], labels = ["Zeit t", "alpha"], labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 10], axesfont = ["HELVETICA", "ROMAN", 8],gridlines=true,axes=boxed);
# Um den Einfluss des Parameters k1 auf die Geschwindigkeit zu zeigen,  plotten wir die Geschwindigkeiten in Abhängigkeit von der Zeit  t  für die Werte  k1 = [ 0.1, 0.5, 1.0] :
> plot([v(t,0.1),v(t,0.5),v(t,1.)],t=0..tW1,title = "\nGeschwindigkeiten eines Körpers bei viskoser Reibung\n", legend = ["k1=0.1", "k1 = 0.5", "k1 = 1.0"], titlefont = ["ARIAL", 15], labels = ["Zeit t", "Geschwindigkeiten"], labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 10], color = [red,blue,green], axesfont = ["HELVETICA", "ROMAN", 8], legendstyle = [font = ["HELVETICA", 9], location = bottom],gridlines=true,axes=boxed);
# Um den Einfluss des Parameters k1 auf die Wurfbahn zu zeigen,  plotten wir die Bahnkurven in Abhängigkeit von der Zeit  t  für die Werte  k1 = [ 0.1, 0.5, 1.0] :
> p1:=plot([x1(t,0.1),x2(t,0.1),t=0..tW1],color =  red):
> p2:=plot([x1(t,0.5),x2(t,0.5),t=0..tW2],color = blue):
> p3:=plot([x1(t,1.0),x2(t,1.0),t=0..tW3],color =green):
# Wir geben abschließend die Wurfbahn aus
> plots[display]([p1,p2,p3],scaling=constrained,title = "\nWurfbahn",titlefont = ["ARIAL", 15], labels = ["x1", "x2"], labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 10], axesfont = ["HELVETICA", "ROMAN", 8],axes=boxed,gridlines=true,color = [red,blue,green]);
# Die Koordinaten der Scheitelpunkte sind:
> x1S:=evalf(x1H(0.1)); x2S:=evalf(x2H(0.1));
> x1S:=evalf(x1H(0.5)); x2S:=evalf(x2H(0.5));
> x1S:=evalf(x1H(1.0)); x2S:=evalf(x2H(1.0));
# und für die größten Wurfweiten erhalten wir
> x1(tW1,0.1);x1(tW2,0.5);x1(tW3,1.0);
> 
;
