% Differenzenverfahren (Aufgabe 19-10):

clear all

a   = 500  ;
b   = 700  ;
EI1 = 6e9  ;
EI2 = 8e10 ;
M   = 800000 ;
c   = 800  ;
q   = 10   ;
k   = 12   ;

L  = a+b    ;           % Laenge
EIB = EI1   ;           % Bezugs-Biegesteifigkeit
nA = 2400   ;           % 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:
iFeder  = round (nA*a/L+ 3.1) ;
iLager  = n - 2 ;

% Belastung
qi(iFeder:iLager) = 0 : q/(iLager-iFeder) : q ;

% Biegesteifigkeit (mue-Werte):
mi(1:iFeder-1) = EI1/EIB ;
mi(iFeder)     = (EI1+EI2)/(2*EIB) ;
mi(iFeder+1:n) = EI2/EIB ;

% Bettungszahl (ki-Werte):
ki(1:iFeder-1) = k   ;
ki(iFeder)     = k/2 + c/h ;
ki(iFeder+1:n) = 0   ;

%Randbedingungen:
A(1:2,:)   = [0 0 0  0 0 1  0 ;
              0 0 0 -1 0 1 0] ;                                         % Einspannung links
A(n-1:n,:) = [    0            0           1      0 0 0 0  ;
              mi(iLager) -2*mi(iLager) mi(iLager) 0 0 0 0] ;     
b(n) = -M*h^2/EIB ;                                                     % Lager und Moment 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

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')   
   
vFeder     = v(iFeder)          % Ausgabe in das Command window
MbLinks    = Mb(3)
MbFeder    = Mb(iFeder)
Mbmax      = max(abs(Mb))
FQLinks    = FQ(3)
Federkraft = c*v(iFeder)