# 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-6:
# 
# Für das obige System sind 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;
                                           b   

                cot(lambda) - a lambda + ------

                                         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);
      1 + b    /  1    \          1        3    /      5\

      ------ + |- - - a| lambda - -- lambda  + O\lambda /

      lambda   \  3    /          45                     
             1 + b    /  1    \          1        3

             ------ + |- - - a| lambda - -- lambda 

             lambda   \  3    /          45        
                                                    (1/2)  

  /                                    (1/2)       \       

1 |        /     2                    \            |       

- \-30 + 6 \225 a  + 150 a + 20 b + 45/      - 90 a/     , 

2                                                          



                                                        (1/2)  

      /                                    (1/2)       \       

    1 |        /     2                    \            |       

  - - \-30 + 6 \225 a  + 150 a + 20 b + 45/      - 90 a/     , 

    2                                                          



                                                      (1/2)  

    /                                    (1/2)       \       

  1 |        /     2                    \            |       

  - \-30 - 6 \225 a  + 150 a + 20 b + 45/      - 90 a/     , 

  2                                                          



                                                        (1/2)

      /                                    (1/2)       \     

    1 |        /     2                    \            |     

  - - \-30 - 6 \225 a  + 150 a + 20 b + 45/      - 90 a/     

    2                                                        
;
> lambda_0:=%[1]; #Startlösung
                                                       (1/2)

     /                                    (1/2)       \     

   1 |        /     2                    \            |     

   - \-30 + 6 \225 a  + 150 a + 20 b + 45/      - 90 a/     

   2                                                        
;
# 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;
                              1.0
                             0.006
                         0.006000000000
                             0.001
                              500.
                              500.
                          0.1666666667
                          1.000000000
;
> 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;
                          1.861736462
;
# und der nummerisch exakte Wert ist
> ew[1]:=fsolve(ewg,lambda=lambda_0);
                          1.814534975
;
# 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);
    Eigenwerte = ,Vector[column](%id = 18446745330023485846)
;
# Wir stellen die ersten n = 3 Eigenfunktionen grafisch dar:
> n:=3; Xn:=seq(sin(ew[i]*xi),i=1..n):
                               3
;
> 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);

> 
;
