% Differenzenverfahren: Gerader Traeger mit Gelenk und Feder (Aufgabe 18-10),
% Berechnung mit der "Sparse-Matrix-Strategie" von MATLAB

clear all

La = 220   ;
Lb = 260   ;
Lc = 300   ;
Ld = 420   ;
EI1 = 2.1e5*2400 ;
EI2 = 2.1e5*800  ;
cF = 300   ;
q0 = 2     ;
F  = 800   ;

L   = La+Lb+Lc+Ld ;     % Laenge
EIB = EI1 ;             % Bezugs-Biegesteifigkeit
nA  = 4800    ;         % Anzahl der Abschnitte
n   = nA + 5  ;         % Anzahl der Gleichungen
h   = L / nA  ;         % Schrittweite

b  = zeros (n,1) ;      % Rechte Seite b
qi = zeros (n,1) ;      % Linienlast
mi = zeros (n,1) ;      % Ii = mi * IB
ki = zeros (n,1) ;      % Bettungszahlen an Stuetzpunkten

% Markante Punkte:
iKraft  = round (nA *  La/L        + 3) ;
iGelenk = round (nA * (La+Lb)/L    + 3) ;
iFeder  = round (nA * (La+Lb+Lc)/L + 3) ;
iRechts = n - 2 ;

% Belastung
qi(iFeder:iRechts) = q0  ;
qi(iKraft)         = F/h ;

% Biegesteifigkeit (mue-Werte):
mi(1:iGelenk-1) = EI1/EIB ;
mi(iGelenk)     = 0       ;
mi(iGelenk+1:n) = EI2/EIB ;

% Feder:
ki(iFeder) = cF/h ;

nRBlinks  = 3 ; % 3 Werte fuer Beschreibung der Einspannung
nRBrechts = 4 ; % 4 Werte fuer Beschreibung der momentenfreien Lagerung
nspare = (n-4)*5+nRBlinks+nRBrechts ;    
AA     = zeros (nspare,1) ;
ipos   = zeros (nspare,1) ;     
jpos   = zeros (nspare,1) ;     

%Einspannung links:
ipos(1:3) = [1  2 2] ;
jpos(1:3) = [3  2 4] ;
AA  (1:3) = [1 -1 1] ;

ii = 4 ;
for i=3:n-2                             % Matrix A:
    ipos(ii:ii+4) = i ;
    jpos(ii:ii+4) = i-2:i+2 ;
    AA  (ii:ii+4) = [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)] ;    
    ii=ii+5 ;
end

% Momentenfreies Lager rechts:
ipos(ii:ii+3) = [n-1 n-1 n-1 n]   ;
jpos(ii:ii+3) = [n-3 n-2 n-1 n-2] ;
AA  (ii:ii+3) = [mi(iRechts) -2*mi(iRechts) mi(iRechts) 1] ;

A = sparse(ipos,jpos,AA,n,n) ; 

% Rechte Seite b:
for i=3:iRechts                          
    b(i) = qi(i)*h^4/(EIB) ;
end

v = A \ b ;    % Loesen des Gleichungssystems (Berechnung der Durchbiegung v)
 
clf;
z = 0 : h : L ;
subplot (3,1,1) ; plot (z , v(3:n-2)) , axis ij , 
                        grid on , title ('Durchbiegung')  % 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')   

format long
AbsenkungKraft  = v(iKraft)   % Ausgabe in das Command window:    
AbsenkungGelenk = v(iGelenk)  % Zahlenwerte der Verschiebungen ...
AbsenkungFeder  = v(iFeder)   % ... an drei speziellen Punkten ...
MbA = Mb(3)                   % ... und das Einspannmoment am Punkt A
