# Kapitel 1: 
# Freie Schwingungen einer rotierenden Saite
# © 2017  Friedrich U. Mathiak, 
# mathiak@mechanik-info.de
# 
> restart: with(LinearAlgebra): with(plots): unprotect(Psi):
# Beispiel 1-11:
# 
# 
# Die freien Schwingungen einer rotierenden Saite werden beschrieben durch die Legendresche Differenzialgleichung
> ode2:=Diff((1-xi^2)*diff(X(xi),xi),xi)+lambda*(lambda+1)*X(xi)=diff((1-xi^2)*diff(X(xi),xi),xi)+lambda*(lambda+1)*X(xi);
# Sie hat die Lösung
> dsolve(rhs(ode2),X(xi)); assign(%);
# Die Legendreschen Funktionen 2. Art Q»(¾) müssen aus physikalischen Gründen ausgeschieden werden, da sie für ¾ ’ 1 über alle Grenzen wachsen. Um die Randbedingung bei ¾ = 0 erfüllen zu können - dort ist die Saite fest eingeklemmt - muss zunächst » = n   ganzzahlig sein. Damit verbleiben als Lösung (auf die Konstante _C1 kommt es hier nicht an) die Legendreschen Polynome:
> X:=unapply(LegendreP(n, xi),n,xi);
> plot([seq(X(i,xi),i=0..5)],xi=0..1,axes=boxed,gridlines=true,title = "\nDie Legendreschen Polynome\n", legend = ["n=0","n=1","n=2","n=3","n=4","n=5"], titlefont = ["ARIAL", 15], labels = ["xi", "LegendreP(n, xi)"], labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 10], linestyle = [solid, dot, dash, dashdot, longdash,spacedash], axesfont = ["HELVETICA", "ROMAN", 8], legendstyle = [font = ["HELVETICA", 9], location = bottom]);
# Um die Randbedingung bei ¾ = 0 erfüllen zu können, dürfen von den Legendreschen Polynomen Pn(¾)  nur die ungeradzahligen verwendet werden. 
> plot([seq(X(2*i-1,xi),i=1..5)],xi=0..1,axes=boxed,gridlines=true,title = "\nDie Legendreschen Polynome (n ungerade)\n", legend = ["n=1","n=3","n=5","n=7","n=9"], titlefont = ["ARIAL", 15], labels = ["xi", "LegendreP(n, xi)"], labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 10], linestyle = [solid, dot, dash, dashdot, longdash], axesfont = ["HELVETICA", "ROMAN", 8], legendstyle = [font = ["HELVETICA", 9], location = bottom]);
# Wir betrachten folgendes Anfangswertproblem: Eine mit der Winkelgeschwindigkeit © um die z-Achse rotierende Saite erhält zum Zeitpunkt t = 0 an der Stelle x = a  einen Hammerschlag senkrecht zur Rotationsebene. Die daraus resultierende Anfangsgeschwindigkeit ¨ wird näherungsweise wie folgt angenommen: 
> Psi:=1/2*v0*(cos(2*Pi*(xi-alpha)/beta)+1);
# Die Fläche, die die Funktion ¨ mit der x-Achse einschließt, ist 
> S:=int(Psi,xi=alpha-beta/2..alpha+beta/2);
> PSI:=unapply(piecewise(xi>=alpha-beta/2 and xi<=alpha+beta/2,Psi),v0,alpha,beta,xi);
# Es gelten die folgenden Systemwerte:
> v0:=10.; l:=1.; a:=l/4; b:= l/5.; Omega:=100.;
> alpha:= a/l; beta:= b/l; 
# Wir stellen die Anfangsgeschwindigkeit grafisch dar:
> plot(PSI(v0,alpha,beta,xi),xi=0..1,axes=boxed,gridlines=true,title = "\nAnfangswert v0\n",titlefont = ["ARIAL", 15], labels = ["¾ = x/", "v0"],labeldirections = ["horizontal", "horizontal"], labelfont = ["HELVETICA", 10]);
# Berechnung der Konstanten Bm, (M: Anzahl der Reihenglieder)
> M:=25: omegam:=Vector(M): Bm:=Vector(M):
> xl:=alpha-beta/2: xr:=alpha+beta/2:  #Integrationsgrenzen
;
> for m to M do
>   omegam[m]:=Omega*sqrt(m*(2*m-1));
>   Bm[m]:=(4*m-1)/omegam[m]*int(Psi*LegendreP(2*m-1,xi),xi=xl..xr);
> end do:
> w:=0: wp:=0:
> for m to M do
>   C := Bm[m]*LegendreP(2*m-1,xi);
>  tau:=omegam[m]*t;
>   w :=w  + C*sin(tau);
>   wp:=wp + C*omegam[m]*cos(tau);
> end do:
> TE:=evalf(2*Pi/omegam[1]);
> para:=xi=0..1,t=0..2*TE,frames=100,thickness=2,gridlines=true,numpoints=150,titlefont = ["ARIAL", 15],labeldirections = ["horizontal", "horizontal"], labelfont = ["HELVETICA", 10],axes=boxed:
> p1:=plots[animate](w,para,title = "\nAuslenkung\n",labels = ["¾ = x/", "w(¾,t)"]): #Das kann etwas dauern!!
;
> p2:=plots[animate](wp,para,title = "\nGeschwindigkeit\n",labels = ["¾ = x/", "v(¾,t)"]): #Das kann etwas dauern!!
;
> plots[display](Matrix(1,2,[p1,p2]));
> 
;
