% Ritzsches Verfahren für die Stabknickung: Aufgabe 33.14

clear all

% 11111111111 ANPASSEN AN DAS AKTUELLE PROBLEM VON HIER ... 111111111111111111111111111
% Parameter:
L1 = 250 ;
L2 = 400 ;
L3 = 150 ;
d1 =   8 ;
d3 =   8 ;
L  = L1 + L2 + L3 ;          
E  = 2.1e5  ;          
c  = 0 ;
cT = 1e10 ;

% Ansatzfunktionen und deren 1. und 2. Ableitungen:
m  = 5 ;                           % Anzahl der Ansatzfuntionen (wenn m < 5 gesetzt ...
P1 = [1/L^3 -2/L^2 1/L   0] ;      % ... gesetzt wird, werden nur die ersten ...
P2 = [1/L^4 -2/L^3 1/L^2 0 0] ;    % ... m Ansatzfunktionen verwendet)
P3 = [1/L^5 -2/L^4 1/L^3 0 0 0] ;
P4 = [1/L^6 -2/L^5 1/L^4 0 0 0 0] ;
P5 = [1/L^7 -2/L^6 1/L^5 0 0 0 0 0] ;
% 11111111111 ... BIS HIER 111111111111111111111111111111111111111111111111111111111111

P1D = polyder (P1);
P2D = polyder (P2);
P3D = polyder (P3);
P4D = polyder (P4);
P5D = polyder (P5);

P1DD = polyder (P1D) ;
P2DD = polyder (P2D) ;
P3DD = polyder (P3D) ;
P4DD = polyder (P4D) ;
P5DD = polyder (P5D) ;

% ------------------------------------------------------------------------------------
% Vorbereitung der numerischen Integration und der Ergebnisauswertung:
n  = 100 ;                % Anzahl der Abschnitte für numerische Integration
dz = L / n ;              % Breite eines Integrationsintervalls
nS = n + 1 ;              % Anzahl der Stützstellen
zS = 0 : dz : L ;         % Koordinaten der Stützpunkte

EIS = zeros (nS , 1) ;
FNQ = zeros (nS , 1) ;

% 22222222222222 ANPASSEN AN DAS AKTUELLE PROBLEM VON HIER ... 222222222222222222222222
kc  = round (n*L1/L + 1) ;               % Angriffspunkt der Feder
kcT = 1 ;                                % Angriffspunkt der Drehfeder
k23 = round (n*(L1+L2)/L + 1) ;          % Übergang Bereich 2 zu 3

d   = zeros (nS , 1) ;
d(1  :kc ) = d1 ;
d(kc :k23) = d1 ;%: (d3-d1)/(k23-kc) : d3 ;
d(k23:nS ) = d3 ;
EIS(1:nS) = E*pi*d.^4/64 ;               % Biegesteifigkeit
FNQ(1:nS) = 1     ;                      % Normalkraft/(Kritische Last)
% 22222222222222 ... BIS HIER 222222222222222222222222222222222222222222222222222222222

% Funktionswerte und 1. und 2. Ableitung der Ansatzfunktionen an den Stützstellen
vS  (:,1) = polyval (P1   , zS)' ;
vS  (:,2) = polyval (P2   , zS)' ;
vS  (:,3) = polyval (P3   , zS)' ;
vS  (:,4) = polyval (P4   , zS)' ;
vS  (:,5) = polyval (P5   , zS)' ;
vdS (:,1) = polyval (P1D  , zS)' ;
vdS (:,2) = polyval (P2D  , zS)' ;
vdS (:,3) = polyval (P3D  , zS)' ;
vdS (:,4) = polyval (P4D  , zS)' ;
vdS (:,5) = polyval (P5D  , zS)' ;
vddS(:,1) = polyval (P1DD , zS)' ;
vddS(:,2) = polyval (P2DD , zS)' ;
vddS(:,3) = polyval (P3DD , zS)' ;
vddS(:,4) = polyval (P4DD , zS)' ;
vddS(:,5) = polyval (P5DD , zS)' ;

subplot(211) ; plot(zS , vS(:,1:m)) ; grid ; title('Ansatzfunktionen') ;

% ------------------------------------------------------------------------------------
% Aufbau des allgemeinen Matrizen-Eigenwertproblems:
K = zeros (m , m) ;
B = zeros (m , m) ;

for ii = 1:m                                          % Schleife über alle Gleichungen
  for jj = 1:m                                        % ii-te Zeile (Koeffizienten kij)
    Summe1 = EIS(1)*vddS(1,ii)*vddS(1,jj) + EIS(nS)*vddS(nS,ii)*vddS(nS,jj) ;
    Summe2 = FNQ(1)*vdS (1,ii)*vdS (1,jj) + FNQ(nS)*vdS (nS,ii)*vdS (nS,jj) ;
    faktor = 4 ;                                      % Numerische Integration ...
    for k = 2:n                                       % ... nach Simpsonscher Regel ... 
      Summe1 = Summe1 + EIS(k) * vddS(k,ii) * vddS (k,jj) * faktor ;
      Summe2 = Summe2 +          vdS (k,ii) * vdS  (k,jj) * faktor ;
      if   (faktor == 4) faktor = 2 ;                 
      else               faktor = 4 ;                 
      end ;
   end
% 44444444444444 ANPASSEN AN DAS AKTUELLE PROBLEM VON HIER ... 4444444444444444444444444
    K(ii,jj) = Summe1*dz/3 + c*vS(kc,ii)*vS(kc,jj) + cT*vdS(kcT,ii)*vdS(kcT,jj) ;
% 44444444444444 ... BIS HIER ... 444444444444444444444444444444444444444444444444444444
    B(ii,jj) = Summe2*dz/3 ;
  end   
end

% Lösen des allgemeinen Eigenwertproblems:
[V D] = eig (K , B) ;
[Fkr , index] = min(diag(D)) ;
FKritisch = Fkr

vSchlange = vS(:,1:m) * V(1:m,index) ;        % Knickfigur zum kleinsten Eigenwert
subplot(212) ; plot(zS , vSchlange) ; grid ; title('Knickfigur') ;