% Differenzenverfahren (Statisch unbestimmter Träger):

clear all

a  = 220 ;
b  = 800 ;
c  = 210 ;
d  = 270 ;
e  = 290 ;
f  = 310 ;
EI = 4.2e9 ;
F1 = 2000 ;
F2 = 1000 ;
q0 = 1 ;

L   = a+b+c+d+e+f ;      % Laenge
EIB = EI    ;            % Bezugs-Biegesteifigkeit
nA  = 2100   ;           % 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:
iB = round (nA*a/L+3) ;
iC = round (nA*(a+b)/L+3) ;
iD = round (nA*(a+b+c)/L+3) ;
iE = round (nA*(a+b+c+d)/L+3) ;
iF = round (nA*(a+b+c+d+e)/L+3) ;
iG = n-2 ;

% Belastung
qi(iB:iC) = 3*q0 : 4*q0/(iC-iB) : 7*q0 ;
qi(iE)    = F1/h ;

% Biegesteifigkeit (mue-Werte):
mi(1:n) = EI/EIB ;

%Randbedingungen:
A(1:2,:)   = [ 0 0 0 0  0 1 0 ;
               0 0 0 1 -2 1 0 ] ;                             % Lager links
A(n-1:n,:) = [ mi(iG-1) -2*mi(iG-1) mi(iG-1)-mi(iG+1) 2*mi(iG+1) -mi(iG+1) 0 0 ;
               mi(iG)   -2*mi(iG)         mi(iG)          0         0      0 0 ] ;     
B(n-1)     = 2*F2*h^3/EIB ;                                   % Kraft F2 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

A(iB,:) = [ 0 0 0 1 0 0 0 ] ;
b(iB)   = 0 ;
A(iC,:) = [ 0 0 0 1 0 0 0 ] ;
b(iC)   = 0 ;
A(iD,:) = [ 0 0 0 1 0 0 0 ] ;
b(iD)   = 0 ;
A(iF,:) = [ 0 0 0 1 0 0 0 ] ;
b(iF)   = 0 ;

tic ;
v = gabamp (A , B) ;                                                    % Berechnung der Durchsenkung
toc ; 

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')   
   
vF1 = v(iE)          % Ausgabe in das Command window
vF2 = v(iG)          
MbB = Mb(iB)
MbC = Mb(iC)
MbD = Mb(iD)
MbE = Mb(iE)
MbF = Mb(iF)
