% Differenzenverfahren, Aufgabe 18-11

clear all

L  = 1200 ;
q0 = 2    ;
EI = 5e9  ;
F1 = 500  ;
F2 = 5000 ;
nA = 600    ;       % Anzahl der Abschnitte (nur dieser Wert muss geaendert werden)
n  = nA + 5 ;       % Anzahl der Gleichungen
h  = L / nA ;

A  = zeros (n,n) ;  % Nullmatrix
b  = zeros (n,1) ;  % Nullvektor
q  = zeros (n,1) ;

A(1:2,1:5)     = [ 0  1 -2  1 0  ; ...
                  -1  2  0 -2 1] ;       % Randbedingungen links: Einzelkraft
b(2) = 2*F1*h^3/EI ;
A(n-1:n,n-4:n) = [ 0  0  1  0 0  ; ...
                   0 -1  0  1 0] ;       % Randbedingungen rechts: Einspannung

iq0  = round (nA/3   + 3.1) ;
i3q0 = round (nA*2/3 + 3.1) ;
q(iq0:i3q0) = q0 : (3*q0 - q0) / (i3q0-iq0) : 3*q0 ;

iF2    = round (nA * 3/4 +3.1) ;         % Kraft F2 "verschmiert"
q(iF2) = F2/h ;
           
for i = 3:n-2                          
    A(i,i-2:i+2)= [1 -4 6 -4 1] ;        % Standardgleichungen
    b(i)        = q(i)*h^4/EI   ;
end

iMitte = round (nA*0.5 + 3.1) ;
A(iMitte,iMitte-2:iMitte+2)= [0 0 1 0 0] ;  % Zwischenstuetze
b(iMitte) = 0 ;

t1 = cputime ;
v = A \ b ;    % Loesen des Gleichungssystems (Berechnung der Durchbiegung v)
ZeitGlSyst = cputime - t1

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

for i = 3:n-2
    Mb(i) = - EI/h^2 * (v(i-1) - 2*v(i) + v(i+1)) ;
    FQ(i) = - EI/(2*h^3) * (-v(i-2) + 2*v(i-1) - 2*v(i+1) + v(i+2)) ;
end

MbMitte  = Mb(iMitte)
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')