# Kapitel 1: 
# Die Saite mit Fesselung des rechten Endes an einem Feder-Masse-System
# © 2017 Friedrich U. Mathiak, 
# mathiak@mechanik-info.de
# 
> restart: with(LinearAlgebra): with(plots):
# Beispiel 1-7:
# 
# Lösung des AWP für die ausgelenkte Saite nach der Kollokationsmethode. Dazu sind vorab die Eigenwerte und Eigenfunktionen zu berechnen. Die Saitenlänge beträgt  = 1,0 m, und die Massendichte ist ñ =7.5 g/cm3.  Die Saite ist mit der Kraft S = 500 N vorgespannt. Die Endmasse M = 0.001 kg  ist an eine lineare Feder mit der Federsteifigkeit k = 500 N/m gefesselt. 
# Wir bestimmen zunächst die Eigenwerte und notieren die Eigenwertgleichung
> ewg:=cot(lambda)-a*lambda+b/lambda;
# Wir beschaffen uns eine Startlösung für den kleinsten Eigenwert. Dazu entwickeln wir die Eigenfunktion in eine Reihe und brechen nach dem 3. Glied ab. Die Auflösung des so entstandenen Polynoms nach » ergibt: 
> series(ewg,lambda,4); p:=convert(%,polynom); solve(p,lambda);
> lambda_0:=%[1]; #Startlösung
;
# Mit den Werten des Beispiels sind:
> l:=1.0; m:=0.006; mu:=m/l; M:=0.001; k:=500.; S:=500.; c:= sqrt(S/mu): a:=M/m; b:=k*l/S;
> plot(ewg,lambda=0..3.*Pi,y=-2..2,gridlines=true,axes=normal);
> nmax:=10: ew:=Vector(nmax):  #Wir berechnen nmax Eigenwerte
;
# Der kleinste Eigenwert ergibt sich mit der obigen Näherungslösung
> lambda_0;
# und der nummerisch exakte Wert ist
> ew[1]:=fsolve(ewg,lambda=lambda_0);
# Wir berechnen die verbleibenden Eigenwerte bis n = nmax
> for i from 2 to nmax do
>   ew[i]:=fsolve(ewg,lambda=ew[i-1]+Pi);
> end do:
> print(`Eigenwerte = `,ew);
# Wir stellen die ersten n = 3 Eigenfunktionen grafisch dar:
> n:=3; Xn:=seq(sin(ew[i]*xi),i=1..n):
> plot([Xn],xi=0..1,title = "\nDie ersten drei Eigenfunktionen\n",titlefont = ["ARIAL", 15], linestyle = [solid, longdash], axes = boxed, axesfont = ["HELVETICA", "ROMAN", 8],legend = ["Grundschwingung", "1. Oberschwingung","2. Oberschwingung"],legendstyle = [font = ["HELVETICA", 9], location = bottom], gridlines=true);
# Lösung des Anfangswertproblems
> h:=0.01;
> Phi:=unapply(h*xi,xi);
> dxi:=1/nmax: xi:=0.: A:=Matrix(nmax): b:=Vector(nmax):
> for i to nmax do
>   xi:=xi + dxi:
>   b[i]:=Phi(xi):
>     for j to nmax do
>       A[i,j]:=sin(ew[j]*xi):
>     end do:
> end do:
> a:=LinearSolve(A,b);
# Näherungslösung für die Anfangsauslenkung ¦:
> xi:='xi':
> Phiq:=add(a[i]*sin(ew[i]*xi/l),i=1..nmax):
> plot(Phiq,xi=0..1,title="\nNäherungslösung für die Anfangsbedingung Phi\n",titlefont = ["ARIAL", 15],gridlines=true);
> xi:='xi':
# Wir führen mit t = T1Ä  =  2*Pi/(omega*`1`);Ä  die  dimensionslose Zeit Ä ein. 
> omega1:=ew[1]*c/l; T1:=2*Pi/omega1; k:=2*Pi/ew[1];
> w:=unapply(add(a[i]*cos(k*ew[i]*tau)*sin(ew[i]*xi),i=1..nmax),xi,tau):
> animate(w(xi,tau),xi=0..1,tau=0..1,title="\nAuslenkung in [m]\n",titlefont = ["ARIAL", 15],frames =100,color=blue,thickness=3,gridlines=true,axes=boxed);
> v:=unapply(T1*diff(w(xi,tau),tau),xi,tau):
> animate(v(xi,tau),xi=0..1,tau=0..1,title="\nGeschwindigkeit in [m/s]\n",titlefont = ["ARIAL", 15],frames =100,color=blue,thickness=3,gridlines=true,axes=boxed);
# Von Interesse ist noch die Verschiebung des rechten Endes als Funktion der Zeit
> plot(w(1,tau),tau=0..5,title="\nAuslenkung des rechten Endes in [m]\n",titlefont = ["ARIAL", 15],gridlines=true,axes=boxed);
> 
;
