# Kapitel 1: 
# Freie Saitenschwingungen mit inhomogener Massenbelegung, 
# Näherungslösung mittels der Störungsrechnung
# © 2017  Friedrich U. Mathiak, 
# mathiak@mechanik-info.de
# 
> restart: with(LinearAlgebra): with(Student[Calculus1]): with(plots): 
# Beispiel 1-10:
# 
# 
# Wir notieren die Partikularlösungen für X0(z):
> XSn[0]:=unapply(cos(z),z); XAn[0]:=unapply(sin(z),z);
# Wir brechen die Reihe nach dem quadratischen Glied ab (N = 2) und beschaffen uns die entsprechenden Partikularintegrale
> N:=2:
> for i to N do
>   XSn[i]:=unapply(int(zeta^2*XSn[i-1](zeta)*sin(z-zeta),zeta=0..z),z):
>   XAn[i]:=unapply(int(zeta^2*XAn[i-1](zeta)*sin(z-zeta),zeta=0..z),z):
> end do:
# Für den symmetrischen und den antimetrischen Anteil der Reihe erhalten wir dann:
> XS:=XSn[0](z): XA:=XAn[0](z):
> for n to N do
>   XS:=XS + epsilon^n*XSn[n](z):
>   XA:=XA + epsilon^n*XAn[n](z):
> end do:
> print(`XS = `,XS); print(`XA = `,XA);
# Mit den Werten des Beispiels:
> l:= 1.; c0:=204.1241452; epsilon:=4*(1-alpha)/kappa0^2; alpha:=1/2; xi:=1/2; z:=kappa0*xi:  
# Nichttriviale Lösungen des homogenen Gleichungssystems (s.h. Buch) existieren nur für DA = 0, also für XS = 0 und XA = 0. Um eine Vorstellung von der Lage der Nullstellen zu bekommen, stellen wir XS und XA grafisch dar. 
> p1:=plot(XS,kappa0=0..35,title = "\nDie Funktion XS\n",titlefont = ["ARIAL",15], axesfont = ["HELVETICA", "ROMAN", 8], gridlines=true,axes=boxed):
> p2:=plot(XA,kappa0=0..35,title = "\nDie Funktion XA\n",titlefont = ["ARIAL",15], axesfont = ["HELVETICA", "ROMAN", 8], gridlines=true,axes=boxed):
> p3:=plot(XS*XA,kappa0=0..35,title = "\nDie Funktion XS*XA\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 XS und XA  im Intervall 0 < º0 <70
> kapp0e:=70;
> EWS:=Roots(XS=0,numeric,1..kapp0e);
> EWA:=Roots(XA=0,numeric,1..kapp0e);
> N:=min(nops(EWS),nops(EWA));
# 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);
# 
