% Differenzenverfahren: Gerader Traeger mit Gelenk und 2 Federn, Loesung
% unter Ausnutzung der Bandstruktur der Koeffizientenmatrix (Aufgabe 18-12)

clear all

La  = 200   ;
Lb  = 700   ;
Lc  = 190   ;
Ld  = 240   ;
Le  = 260   ;
Lf  = 280   ;
EI1 =  6e9  ;
EI2 = 12e9  ;
EI3 =  1e9  ;
c1 =  100   ;
c2 =   50   ;
q1 =    3   ;
F  = 2000   ;
M  = 200000 ;

L  = La+Lb+Lc+Ld+Le+Lf ; % Länge
EIB = EI1 ;              % Bezugs-Biegesteifigkeit
nA  = 1870    ;          % Anzahl der Abschnitte
n   = nA + 5  ;          % Anzahl der Gleichungen
h   = L / nA  ;          % Schrittweite

A  = zeros (n,7) ;      % Rechteckmatrix fuer Aufnahme der Bandmatrix
b  = zeros (n,1) ;      % Nullvektor ("rechte Seite")
qi = zeros (n,1) ;      % Linienlastintensitaeten an Stuetzpunkten
mi = zeros (n,1) ;      % EIi = mi * EIB
ki = zeros (n,1) ;      % Bettungszahlen an Stuetzpunkten

% Markante Punkte:
iq1     = round (nA *  La/L              + 3) ;
iG      = round (nA * (La+Lb)/L          + 3) ;
ic1     = round (nA * (La+Lb+Lc)/L       + 3) ;
iF      = round (nA * (La+Lb+Lc+Ld)/L    + 3) ;
iLager  = round (nA * (La+Lb+Lc+Ld+Le)/L + 3) ;
ic2     = n - 2 ;

% Belastung
qi(iq1:iG) = q1 : -q1/(iG-iq1) : 0 ;
qi(iF)     = F/h ;

% Biegesteifigkeit (mue-Werte):
mi(1:iq1-1)    = EI1/EIB ;
mi(iq1)        = (EI1+EI2)/(2*EIB) ;
mi(iq1+1:iF-1) = EI2/EIB ;
mi(iG)         = 0       ;
mi(iF)         = (EI2+EI3)/(2*EIB) ;
mi(iF+1:n)     = EI3/EIB ;

% Feder c1:
ki(ic1) = c1/h ;

%Randbedingungen:
A(1:2,:)           = [0 0  0   0      0       1   0  ;
                      0 0  0 mi(3) -2*mi(3) mi(3) 0] ;
b(2)                 = M*h^2/EIB ;                           % Lager und Moment links
A(n-1:n,:) = [mi(n-3) -2*mi(n-3) mi(n-3)-mi(n-1)+2*c2*h^3/EIB 2*mi(n-1) -mi(n-1)  0 0  ;
              mi(ic2) -2*mi(ic2)              mi(ic2)             0         0     0 0] ;
                        % Feder rechts     

for i = 3:n-2    % Matrix A:                          
    A(i,:)= [0 mi(i-1) -2*(mi(i-1)+mi(i)) mi(i-1)+4*mi(i)+mi(i+1)+ki(i)*h^4/(EIB) ...
                       -2*(mi(i)+mi(i+1)) mi(i+1) 0] ;    
    b(i) = qi(i)*h^4/(EIB) ;                          % Standardgleichungen
end

%Lager:
A(iLager,:)= [0 0 0 1 0 0 0] ;
b(iLager)  = 0 ;

t1 = cputime ;
v = gabamp (A , b) ;                                  % Berechnung der Durchsenkung
ZeitGlSyst = cputime - t1 

clf;
z = 0 : h : L ;
subplot (3,1,1) ; plot (z , v(3:n-2)) , axis ij , grid on , 
                        title ('Durchbiegung')  % Graphische Ausgabe der Biegelinie

for i=3:n-2                       
    Mb(i) = -mi(i)*EIB * (v(i-1)-2*v(i)+v(i+1))/h^2 ;                 % Biegemoment
    FQ(i) = 0.5*EIB * (mi(i-1)*v(i-2)-2*mi(i-1)*v(i-1)+(mi(i-1)-mi(i+1)) ...
                        *v(i)+2*mi(i+1)*v(i+1)-mi(i+1)*v(i+2))/h^3 ;  % Querkraft
end
    
subplot (3,1,2) ; plot (z , Mb(3:n-2)) , grid on , title ('Biegemoment') 
subplot (3,1,3) ; plot (z , FQ(3:n-2)) , grid on , title ('Querkraft')   
   
AbsenkungGelenk  = v(iG) 
AbsenkungKraft   = v(iF)         % Ausgabe in das Command window:    
AbsenkungFeder1  = v(ic1)        % Zahlenwerte der Verschiebungen ...
AbsenkungFeder2  = v(ic2)        % an vier speziellen Punkten ...
Mbmax            = max(abs(Mb))  % ... und das absolut grösste Biegemoment
