clear
clc

%Biegungsgleichung : FDM
% H = EIy'      : [-1 +1]/dx
% M = EIy''     : [ 1 -2 1]/dx^2
% V = EIy'''    : [-1 +3 -3 1]/dx^3
% w = EIy''''   : [ 1 -4 6 -4 1]/dx^4


%% Eingangswerte
%Eingangswerte Balken
L = 1;          %Länge des Balkens  [m]
B = 0.1;        %Breite des Balkens [m]
H = 0.1;        %Höhe des Balkens   [m]
I = B*(H^3)/12;  %Trägheitsmoment 2.Grades   [m^4]
E = 210e9;      %E-Modul Stahl              [N/m^2]
EI = E*I;       %Biegesteifigkeit           [N*m^2]
rhoA = 3;
%Belastung 
w0 = 10e3;        %Streckenlast       [N/m]
%Diskretisierung 
ne = 16;        %Anzahl der Elemente
nx = ne + 1;    %Anzahl der Knoten
dx = L/ne;      %Länge der Elemente
X = 0:dx:L;
%Eigenschwingung
nA = 6;
rhoA = 3;
h = L/nA;
n = nA-1;

%% Matrizen
f = zeros(nx,1);            %Lastvektor
K = zeros(nx,nx);           %Steifigkeitsmatrix
y = zeros(nx,1);            %Biegelinie
P = zeros(nx,1);            %Belastung


% Differentialgleichung (GDE)
% GDE: EIy'''' = -q
% FDM: EI[y(i-2) -4y(i-1) +6y(i) -4y(i+1) +y(i+2)]/dx^4 = -q
%    : EI[y(i-2) -4y(i-1) +6y(i) -4y(i+1) + y(i+2)] = -q*dx^4
%    : EI[ 1 -4 6 -4 1 ] in K
%
% [K][y] = [f]          %Elementsteifigkeitsbeziehung

%% Lastvektor
f = ones(nx,1)*-w0*dx^4/EI;     %Lastvektor

%% Steifigkeitsmatrix
for i = 3:nx-2
    K(i,i-2) =  1;
    K(i,i-1) = -4;
    K(i,i  ) =  6;
    K(i,i+1) = -4;
    K(i,i+2) =  1;
end


%Randbedingungen
%% Auflager links         
K(1,1) =  1;           % left y_1 = 0
f(1) = 0;

%Moment links
% M_1 = EIy'' = 0
% FDM: EI[y(1) -2y(2) +y(3)]/dx^2 = 0
%    : [ 1 -2  1 ] in K
K(2,1) =  1;           % links M_1 = 0
K(2,2) = -2;   
K(2,3) =  1;  

%% Auflager rechts 
K(nx  ,nx  ) =  1;     % rechts y_n = 0
f(nx) = 0;

%Moment rechts
% M_n = EIy'' = 0
% FDM: EI[y(n-2) -2y(n-1) +y(n)]/dx^2 = 0
%    : [ 1 -2  1 ] in K
K(nx-1,nx-2) =  1;     % rechts M_n = 0
K(nx-1,nx-1) = -2; 
K(nx-1,nx  ) =  1;

%% Inverse und Lösung
y = K\f;                %Verschiebung
Ki=inv(K);              %Inverse K = Nachgiebigkeit

%%K anpassen
Keig = K(3:end-2,3:end-2); % Auflager raus
 
Keig=Keig(1:2:end,1:2:end); % verdrehung raus


%%Eigenwertzerlegung K
[V,D] = eig(Keig);

ds=L/(size(V,2)-1);
xx=0:ds:L;

Npp = ceil(sqrt(size(V,2)));

figure
for ee=1:size(V,2)
     
  subplot(Npp,Npp,ee)
  plot(xx(1:end) , V(1:end,ee) )
  title(['f=' num2str( sqrt(D(ee,ee))/2/pi )]);
       
 end