% Differenzenverfahren, Aufgabe 18-11, Loesung
% unter Ausnutzung der Bandstruktur der Koeffizientenmatrix

clear all

a  = 400 ;
b  = 200 ;
c  = 100 ;
d  = 300 ;
L  = a+2*b+c+d ;
q0 = 2    ;
EI = 5e9  ;
F1 = 500  ;
F2 = 5000 ;

nA   = 4800   ;     % Anzahl der Abschnitte (nur dieser Wert muss geaendert werden)
n    = nA + 5 ;     % Anzahl der Gleichungen
h    = L / nA ;
EIB  = EI ;         % Bezugs-Biegesteifigkeit

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:
iq0    = round (nA*a/L         + 3) ;
i3q0   = round (nA*(a+2*b)/L   + 3) ;
iF2    = round (nA*(a+2*b+c)/L + 3) ;  
iLager = round (nA*(a+b)/L     + 3) ;

% Belastung:
qi(iq0:i3q0) = q0 : (3*q0 - q0) / (i3q0-iq0) : 3*q0 ;
qi(iF2) = F2/h ;     % Kraft F2 "verschmiert"

% Biegesteifigkeit (mue-Werte):
mi(1:n) = EI/EIB ;
           
A(1:2,:) = [ 0 0  0       0       mi(3)    -2*mi(3)  mi(3) ;
             0 0 mi(2) -2*mi(2) mi(2)-mi(4) 2*mi(4) -mi(4) ] ; 
B(2) = -2*F1*h^3/EI ;                       % Randbedingungen links: Einzelkraft
A(n-1:n,:) = [ 0 0 1 0 0 0 0 ;
              -1 0 1 0 0 0 0 ] ;            % Randbedingungen rechts: Einspannung

for i = 3:n-2                          
    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/EI   ;              % Standardgleichungen
end

A(iLager,:) = [0 0 0 1 0 0 0] ;                 % Zwischenstuetze
B(iLager)   = 0 ;

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

format long
vF1 = v(3)
vF2 = v(iF2) 

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

MbLager  = Mb(iLager)
MbRechts = Mb(n-2)

z = 0 : h : L ;
subplot (3 , 1 , 1) ; plot (z , v (3:n-2)) , axis ij , grid on , title ('Verschiebung')
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')