# Kapitel 1: 
# Saite mit diskreten Massen
# © 2017  Friedrich U. Mathiak, 
# mathiak@mechanik-info.de
# 
> restart:with(LinearAlgebra): with(plots): with(plottools):
> Proc_Calc_01:=proc(Daten::list)
> #----------------------------------------------------------------------------
> #Eingabe:
> #   Daten[1] = [m,S,l,N]
> #               m: Teilmasse         [kg]
> #               S: Vorspannkraft      [N]
> #               l: Saitenlänge        [m]
> #               N: Anzahl Teilmassen
> #   Daten[2] = [y01,y02,...,y0N] = y0: Anfangsauslenkungen      in [m]
> #   Daten[3] = [v01,v02,...,v0N] = v0: Anfangsgeschwindigkeiten in [m/s]
> #Ausgabe: (Sämtliche global vereinbarten Größen)
> #           omega: Eigenkreisfrequenzen [rad/s], Vektor(N)
> #               f: Frequenzen [Hz], Vektor(N)
> #               A: Amplituden, Vektor(N)
> #           alpha: Phasenverschiebungswinkel
> #----------------------------------------------------------------------------
> local m,S,l,N,dl,y0,v0,ID,i,j,ej,ejT,sina,cosa;
> global omega,f,A,alpha,K,EW,Phi,y,v;
> description "Die beidseitig eingeklemmte Saite mit diskret verteilter Masse";
> S:=Daten[1,2];l:=Daten[1,3]; N:=Daten[1,4]; dl:=l/(N+1); m:=Daten[1,1]/(N+1);
> if N < 1 then
>   print(`Anzahl N der Teilmassen >= 1 wählen`);
>   return; 
> end if:
> y0:=Vector(Daten[2]); v0:=Vector(Daten[3]);
> K :=Matrix(N,N,0);    ID:=LinearAlgebra[IdentityMatrix](N);
> omega:=Vector(N); f:=Vector(N); A:=Vector(N); alpha:=Vector(N);
> #Aufstellen der Systemmatrix
> for i to N-1 do
>   K[i  ,  i]:= 2.;
>   K[i  ,i+1]:=-1.;
>   K[i+1,  i]:=-1.;
> end do;
> K[N,N]:=2.;
> #Berechnung der Eigenwerte und Eigenvektoren
> EW,Phi:=LinearAlgebra[Eigenvectors](K); 
> EW:=Re(EW);Phi:=Re(Phi);
> #Berechnung der Eigenkreisfrequenzen, Kreisfrequenzen und Amplituden Aj
> for j to N do
>   omega[j]:=sqrt(EW[j]*S/m/dl);
>   f[j]:=omega[j]/2/Pi; 
>   ej  :=LinearAlgebra[Column](Phi,j);
>   ejT :=LinearAlgebra[Transpose](ej);
>   A[j]:=sqrt((ejT.y0)^2 + (ejT.v0/omega[j])^2)/(ejT.ej);
>     if abs(A[j]) > 10^(-7) then
>       sina:=ejT.v0/omega[j];
>       cosa:=ejT.y0;
>       alpha[j]:=arctan(sina,cosa);
>     end if;
> end do;
> #Lösung
> y:=0:
> for i to N do
>   y:= y + A[i]*LinearAlgebra[Column](Phi,i)*cos(omega[i]*t-alpha[i]);
> end do;
>   v:=map(diff,y,t);
> end proc:
> 
;
# Beispiel 1-1:
# 
# Auf der in der obigen Abbildung skizzierten elastischen Saite befinden sich im Abstand d drei symmetrisch verteilte frei bewegliche Teilmassen m. Die Saite ist mit der Spannkraft S = 55 N vorgespannt und an den Enden eingeklemmt, womit w0 = w4 = 0 zu fordern sind. Es sind folgende Teilaufgaben zu lösen:
# 1. Berechnung der Eigenwerte und Eigenvektoren. 
# 2. Berechnung des Bewegungsvorgangs, wenn die Massen ohne Anfangsgeschwindigkeit aus der Ruhelage mit w10 = w20 = w30 = 1 cm ausgelenkt und dann sich selbst überlassen werden. 
# 3. Stellen sie zur automatisierten Auswertung der Bewegungsgleichung eine Maple-Prozedur zur Verfügung, und animieren Sie den Bewegungsvorgang. 
> unprotect(D);
# Mit den Werten des Beispiels sind
> M:= 0.4594579258e-3;   #Gesamtmasse in [kg] 
;
> l:=0.65;               #Saitenlänge in [m]
;
> S:=55.;                #Spannkraft  in [N]
;
> N:=3;                  #Anzahl der Teilmassen
;
> d:=l/(N+1);            #Abstand der Teilmassen in [m]
;
> m:= M/(N+1);           #Teilmasse in [kg]
;
> sys:=[M,S,l,N]: 
> y0:=[0.01,0.01,0.01]:  #Anfangsauslenkungen      in   [m]
;
> v0:=[0,0,0]:           #Anfangsgeschwindigkeiten in [m/s]
;
> eingabe1:=[sys,y0,v0];
> Proc_Calc_01(eingabe1):
> EW; Phi; y; v;
# 
> Proc_Calc_02:=proc(Daten)
> #----------------------------------------------------------------------------
> #Eingabe: Daten[1]: Länge der Saite l in [m]
> #         Daten[2]: Anzahl der Teilmassen N
> #         Daten[3]: Auslenkung y in [m]
> #         Daten[4]: Simulationszeit te in [s]
> #         Daten[5]: Anzahl der Teilbilder (frames)
> #Ausgabe:
> #         Animation des Bewegungsvorganges AN 
> #----------------------------------------------------------------------------
> local N,d,par1,par2,P,G,g,i,p,la,le;
> global AN;
> description "Animation der beidseitig eingeklemmten Saite mit diskret verteilter Masse";
> N:=Daten[2]; d:=Daten[1]/(N+1);
> par1:= symbol = solidcircle, symbolsize = 20, color = black:
> par2:= t = 0 .. Daten[4], frames = Daten[5]:
> P:={}: G:={}:
> for i to N do
>   p:=animate(pointplot,[[i*d, y[i]], par1], par2):
>   P:=P union {p};
> end do:
> for i to N-1 do
>   g:=animate(polygonplot,[[[i*d,y[i]],[(i+1)*d,y[i+1]]],thickness = 2, color = red],par2):
>   G:=G union {g};
> end do:
>   la:=animate(polygonplot,[[[  0,   0],[      d,y[1]]],thickness = 2,color = red],par2):
>   le:=animate(polygonplot,[[[N*d,y[N]],[(N+1)*d,   0]],thickness = 2,color = red],par2):
>   G :=G union {la,le}:
> AN:=display(P,G,gridlines,axes=boxed,title = "\nDie Saite mit diskret verteilter Masse\n", titlefont = ["ARIAL", 15], labels = ["Länge der Saite in [m]", "Auslenkung in [m]"], labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 10], axesfont = ["HELVETICA", "ROMAN", 8]);
> end proc:
> te:=0.05; nrahmen:=150:
> eingabe2:=[l,N,y,te,nrahmen]:
> Proc_Calc_02(eingabe2); #Animation des Bewegungsvorgangs
;
# Wir erhöhen die Anzahl der Teilmassen  auf  N = 7  und verändern die Anfangsbedingungen:
> N:=7:
> sys:=[M,S,l,N]: 
> y0:=[0,0,0.01,0.01,0,0,0]:   #Anfangsauslenkungen      in   [m]
;
> v0:=[1.,0,0,0,0,0,0]:        #Anfangsgeschwindigkeiten in [m/s]
;
> eingabe1:=[sys,y0,v0];
> Proc_Calc_01(eingabe1):      #Berechnung der Auslenkungen
;
> eingabe2:=[l,N,y,te,nrahmen]:
> Proc_Calc_02(eingabe2);      #Animation des Bewegungsvorgangs
;
> 
;
