% Differenzenverfahren, Beispiel auf den Seiten 260 bis 262,
% Ausnutzung der Bandstruktur der Koeffizientenmatrix

clear all

L  = 1 ;
q0 = 1 ;
EI = 1 ;
nA = 1000   ;       % Anzahl der Abschnitte (nur dieser Wert muss geaendert werden)
n  = nA + 5 ;       % Anzahl der Gleichungen
h  = L / nA ;

hdsp = 4 ;          % Hauptdiagonalelemente stehen in der 4. Spalte der Bandmatrix A

A  = zeros (n,7) ;  % Nullmatrix fuer Aufnahme der Bandmatrix
b  = zeros (n,1) ;  % Nullvektor

A(1:2,:)   = [ 0  0  0  0  0  1 0 ;
               0  0  0 -1  0  1 0 ] ;     % Randbedingungen links: Einspannung
A(n-1:n,:) = [-1  2  0 -2  1  0 0 ;
               1 -2  1  0  0  0 0 ] ;     % Randbedingungen rechts: Freier Rand

for i = 3:n-2                          
    A(i,:) = [0 1 -4 6 -4 1 0] ;            % Standardgleichungen
    b(i)   = q0*h^4/EI         ;
end

% Lösen des Gleichungssystems 
tic ;
v = gabamp (A , b , hdsp) ;               % Berechnung der Durchsenkung
toc

format long
vEnd    = v(n-2)
phiEnd  = (-v(nA+2) + v(nA+4)) / (2*h)
Mblinks = -EI/h^2 * (v(2) - 2*v(3) + v(4))

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

z = 0 : h : L ;
subplot (2 , 1 , 1) ; plot (z , v (3:n-2)) , axis ij , title ('Verschiebung')
subplot (2 , 1 , 2) ; plot (z , Mb(3:n-2)) , title ('Biegemoment')