# Kapitel 1: 
# Freie Saitenschwingungen mit inhomogener Massenbelegung
# © 2017  Friedrich U. Mathiak, 
# mathiak@mechanik-info.de
# 
> restart: with(LinearAlgebra): with(Student[Calculus1]): with(plots): 
# Beispiel 1-9:
# 
# Für eine Saite mit quadratisch verteilter Massenbelegung ¼ sind die Eigenwerte, Eigenkreisfrequenzen und die Eigenfunktionen zu bestimmen. 
> mu:=unapply(mu0*(4*xi*(xi-1)*(alpha-1)+alpha),mu0,alpha,xi);
> plot(mu(1,1/2,xi),xi=0..1, labels=[xi, mu],title = "\nMassenbelegung der Saite\n",titlefont = ["ARIAL",15], gridlines=true,axes=boxed,filled = [color = "Blue", transparency = .5]);
> ode2:=simplify(diff(X(xi),xi,xi)+kappa0^2*mu(mu0,alpha,xi)/mu0*X(xi));
> dsolve(ode2,X(xi)); assign(%);
> X:=unapply(X(xi),kappa0,alpha,xi,_C1,_C2);
# Wir passen die Verschiebungsfunktion X an die Randbedingung der beidseitig eingeklemmten Saite an:
> gl1:=X(kappa0,alpha,0,_C1,_C2); gl2:=X(kappa0,alpha,1,_C1,_C2);
# und lösen das homogene Gleichungssystem nach den Unbekannten _C1 und _C2 auf.
> gl :=[gl1,gl2]: unb:=[_C1,_C2]:
> A,b:=GenerateMatrix(gl,unb):
> DA :=Determinant(A);
# Die Determinante ergibt sich als das Produkt der beiden Funktionen
> H1:=unapply(op(2,DA),kappa0,alpha); H2:=unapply(op(3,DA),kappa0,alpha);
# Mit den Werten des Beispiels können wir die Aufgabenstellung weiter konkretisieren:
> l:=1.0; mu0:=0.012; alpha:=1/2; S:=500.; c0:= sqrt(S/mu0);
> X(kappa0,alpha,xi,_C1,_C2);
> X(kappa0,alpha,0,_C1,_C2); X(kappa0,alpha,1,_C1,_C2);
> H1(kappa0,alpha); H2(kappa0,alpha);
# Damit erhalten wir:
> DA;
> H1(kappa0,alpha);H2(kappa0,alpha);
# Nichttriviale Lösungen des homogenen Gleichungssystems existieren nur für DA = 0, also für H1 = 0 und H2 = 0. Um eine Vorstellung von der Lage der Nullstellen zu bekommen, stellen wir H1 und H2 und das Produkt H1*H2 grafisch dar. 
> p1:=plot(H1(kappa0,alpha),kappa0=0..35,title = "\nDie Funktion H1\n",titlefont = ["ARIAL",15], axesfont = ["HELVETICA", "ROMAN", 8], gridlines=true,axes=boxed):
> p2:=plot(H2(kappa0,alpha),kappa0=0..35,title = "\nDie Funktion H2\n",titlefont = ["ARIAL",15], axesfont = ["HELVETICA", "ROMAN", 8],gridlines=true,axes=boxed):
> p3:=plot(H1(kappa0,alpha)*H2(kappa0,alpha),kappa0=0..35,title = "\nDie Funktion H1*H2\n",titlefont = ["ARIAL",15], axesfont = ["HELVETICA", "ROMAN", 8], gridlines=true,axes=boxed):
> display(Matrix(1,3,[p1,p2,p3]));
# Wir ermitteln die  Nullstellen (Eigenwerte) der Funktionen H1 und H2  im Intervall 0 < º0 <70
> EWS:=Roots(H1(kappa0,alpha)=0,numeric,0..70);
> EWA:=Roots(H2(kappa0,alpha)=0,numeric,0..70);
> N:=min(nops(EWS),nops(EWA));
# Mit den nummerisch ermittelten Eigenwerten ergeben sich folgende Eigenfunktionen:
> EFS:=unapply(X(kappa0,alpha,xi,1,0),kappa0,xi); EFA:=unapply(X(kappa0,alpha,xi,0,1),kappa0,xi); 
# Wir stellen jeweils die ersten drei Eigenfunktionen mit den Eigenwerten EWS und EWA grafisch dar:
> para1:=xi=0..1,titlefont = ["ARIAL",15],axesfont = ["HELVETICA", "ROMAN", 8],legend = ["kappa01","kappa02","kappa03"],legendstyle = [font = ["HELVETICA", 9],location=bottom],gridlines=true,axes=boxed:
> p1:=plot([seq(EFS(EWS[i],xi),i=1..3)],para1,title = "\nDie symmetrischen Eigenfunktion EFS\n"):
> p2:=plot([seq(EFA(EWA[i],xi),i=1..3)],para1,title = "\nDie antimetrischen Eigenfunktion EFA\n"):
> display(Matrix(1,2,[p1,p2]));
# Die Eigenfunktionen EFS und EFA symmetrisch bzw. antimetrisch bezüglich ¾ = 1/2 und erfüllen die geforderten Randbedingungen  bei ¾ = 0 und ¾ = 1.
# Mit den Eigenwerten º0 liegen auch die Eigenkreisfrequenzen É = º0 c0 /  fest. 
> OMS:=c0*EWS/l: OMA:=c0*EWA/l:
> print(`OMS = `,OMS);
> print(`OMA = `,OMA);
> 
;
# Im nächsten Schritt soll der Bewegungsvorgang beschrieben werden. Dazu benötigen wir den zeitabhängigen Lösungsteil aus dem Bernoullischen Produktansatz.
> h:=0.01: a:=0.3:
> PHI:=unapply(piecewise(xi>=0 and xi<=a/l,h/a*xi*l,xi>=a/l and xi<=1, h*(1-xi)/(1-a/l)),xi):
> plot(PHI(xi),xi=0..1,title = "\nAnfangsauslenkung Phi der gezupften Saite\n",titlefont = ["ARIAL",15],axesfont = ["HELVETICA", "ROMAN", 8],gridlines=true,axes=boxed);
> X11:=Matrix(N): X12:=Matrix(N): X21:=Matrix(N): X22:=Matrix(N): Phi:=Vector(2*N):
# Festlegung der Stützpunkte und Funktionswerte der Anfangsauslenkung
> XI:=seq(l/(2*N+1)*i,i=1..2*N): for i to 2*N do Phi[i]:=PHI(XI[i]) end do:
# Aufbau und Lösung der Kollokationsgleichungen
> for i to N do 
>   for j to N do
>     X11[i,j]:=evalf(EFS(EWS[j],XI[i  ]));
>     X12[i,j]:=evalf(EFA(EWA[j],XI[i  ]));
>     X21[i,j]:=evalf(EFS(EWS[j],XI[i+N]));
>     X22[i,j]:=evalf(EFA(EWA[j],XI[i+N]));
>   end do:
> end do:
> An:=LinearSolve(Matrix([[X11,X12],[X21,X22]]),Phi):
# Wir führen zur Animation des Bewegungsvorganges die dimensionslose Zeit Ä ein.
> T1:=2*Pi*l/EWS[1]/c0: t:=tau*T1:
> w:=0.: xi:='xi':
# Das führt auf die Verschiebung w(¾,Ä)
> for i to N do
>   w:= w + evalf(EFS(EWS[i],xi)*An[i]*cos(OMS[i]*t) + EFA(EWA[i],xi)*An[i+N]*cos(OMA[i]*t)):
> end do:
# Durch Ableitung nach der Zeit t erhalten wir die Geschwindigkeit v(¾,Ä)
> v:=diff(w/T1,tau):
# Wir animieren den Bewegungsvorgang:
> para2:=xi=0..1,tau=0..1,titlefont=["ARIAL", 15],frames =100,color=blue,thickness=3,gridlines=true,numpoints=100,axes=boxed:
> an1:=animate(w,para2,title="\nAuslenkung in [m]\n"):
> an2:=animate(v,para2,title="\nGeschwindigkeit in [m/s]\n"):
> display(Matrix(1,2,[an1,an2]));
> 
;
