% Autor: S³awomir Milewski, PK WIL L5, slawek@L5.pk.edu.pl
clear all
close all
% TYPY PODPARCIA BELKI
% 1 - podpora swobodna (blokada ugiêcia po 'y'), 
% 2 - utwierdzenie (blokada ugiêcia po 'y' i obrotu po 'z')
% 3 - utwierdzenie z przesuwem po osi 'y' (blokada obrotu po 'z')
% 4 - koniec swobodny (lub podparcie przegubowe przesuwne po osi 'y')

typ1 = 1; % typ podpory na lewym koñcu belki
typ2 = 1; % typ podpory na prawym koñcu belki
winkler = 0; % w³¹czenie / wy³¹czenie z obliczeñ pod³o¿a sprê¿ystego Winklera
b = 1; % szerokoœæ przekroju - TYLKO DLA WINKLERA
c = 1; % modu³ podatnoœci pod³o¿a - TYLKO DLA WINKLERA

% komplet wzorów ró¿nicowych IV rzêdu dla pochodnych "0",I,II,III i IV -
% wyprowadzenie z szeregu Taylora
P = [1 -2 0.5*2^2 -(1/6)*2^3 (1/24)*2^4;
    1 -1 0.5 -(1/6) (1/24);
    1 0 0 0 0;
    1 1 0.5 (1/6) (1/24);
    1 2 0.5*2^2 (1/6)*2^3 (1/24)*2^4];
W = P'^-1 * eye(5,5);
format rat
W
format short g

n = 7; % liczba wêz³ów
q = 3; % obci¹¿enie ci¹g³e (na ca³ej d³ugoœci belki)
L = 3; % d³ugoœæ belki
EJ = 2e4; % sztywnoœæ na zginanie
P1 = 0; % si³a skupiona na lewym koñcu belki
P2 = 0; % si³a skupiona na prawym koñcu belki
M1 = 0; % moment skupiony na lewym koñcu belki
M2 = 0; % moment skupiony na prawym koñcu belki

h = L/(n-1); % modu³ siatki

A = zeros(n+4,n+4); % macierz wspó³czynników
B = zeros(n+4,1); % wektor prawej strony

% kollokacja we wêz³ach obszaru
for i=1:n
    A(i+2,i+2) = 6/(h^4) + (b*c/EJ) * winkler;
    A(i+2,i+1) = -4/(h^4);
    A(i+2,i+3) = -4/(h^4);
    A(i+2,i) = 1/(h^4);
    A(i+2,i+4) = 1/(h^4);
    B(i+2) = q/EJ;
end
% kollokacja we wêŸle brzegowym lewym
if (typ1==1)
    A(1,3) = 1;
    A(2,1) = -(1/12)/(h^2);
    A(2,2) = (4/3)/(h^2);
    A(2,3) = -(5/2)/(h^2);
    A(2,4) = (4/3)/(h^2);
    A(2,5) = -(1/12)/(h^2);
    B(2) = -M1/EJ;
end
if (typ1==2)
    A(1,3) = 1;
    A(2,1) = (1/12)/h;
    A(2,2) = (-2/3)/h;
    A(2,3) = 0;
    A(2,4) = (2/3)/h;
    A(2,5) = -(1/12)/h;
end
if (typ1==3)
    A(1,1) = (1/12)/h;
    A(1,2) = (-2/3)/h;
    A(1,3) = 0;
    A(1,4) = (2/3)/h;
    A(1,5) = -(1/12)/h;
    
    A(2,1) = -(1/2)/(h^3);
    A(2,2) = 1/(h^3);
    A(2,3) = 0;
    A(2,4) = -1/(h^3);
    A(2,5) = (1/2)/(h^3);
    B(2) = P1/EJ;
end
if (typ1==4)
    A(1,1) = -(1/12)/(h^2);
    A(1,2) = (4/3)/(h^2);
    A(1,3) = -(5/2)/(h^2);
    A(1,4) = (4/3)/(h^2);
    A(1,5) = -(1/12)/(h^2);
    B(1) = -M1/EJ;
    
    A(2,1) = -(1/2)/(h^3);
    A(2,2) = 1/(h^3);
    A(2,3) = 0;
    A(2,4) = -1/(h^3);
    A(2,5) = (1/2)/(h^3);
    B(2) = P1/EJ;
end
% kollokacja we wêŸle brzegowym prawym
if (typ2==1)
    A(n+4,n+2) = 1;
    A(n+3,n+4) = -(1/12)/(h^2);
    A(n+3,n+3) = (4/3)/(h^2);
    A(n+3,n+2) = -(5/2)/(h^2);
    A(n+3,n+1) = (4/3)/(h^2);
    A(n+3,n) = -(1/12)/(h^2);
    B(n+3) = -M2/EJ;
end
if (typ2==2)
    A(n+4,n+2) = 1;
    A(n+3,n+4) = -(1/12)/h;
    A(n+3,n+3) = (2/3)/h;
    A(n+3,n+2) = -0;
    A(n+3,n+1) = -(2/3)/h;
    A(n+3,n) = (1/12)/h;
end
if (typ2==3)
    A(n+4,n+4) = -(1/12)/h;
    A(n+4,n+3) = (2/3)/h;
    A(n+4,n+2) = 0;
    A(n+4,n+1) = -(2/3)/h;
    A(n+4,n) = (1/12)/h;

    A(n+3,n+4) = -(1/2)/(h^3);
    A(n+3,n+3) = 1/(h^3);
    A(n+3,n+2) = 0;
    A(n+3,n+1) = -1/(h^3);
    A(n+3,n) = (1/2)/(h^3);
    B(n+3) = P2/EJ;
end

if (typ2==4)
    A(n+4,n+4) = -(1/12)/(h^2);
    A(n+4,n+3) = (4/3)/(h^2);
    A(n+4,n+2) = -(5/2)/(h^2);
    A(n+4,n+1) = (4/3)/(h^2);
    A(n+4,n) = -(1/12)/(h^2);
    B(n+4) = -M2/EJ;
    
    A(n+3,n+4) = -(1/2)/(h^3);
    A(n+3,n+3) = 1/(h^3);
    A(n+3,n+2) = 0;
    A(n+3,n+1) = -1/(h^3);
    A(n+3,n) = (1/2)/(h^3);
    B(n+3) = P2/EJ;
end
A
B
d = det(A);
eps = 10^-16;
if (abs(d)<eps)
    disp('KONSTRUKCJA GEOMETRYCZNIE ZMIENNA!!! - ZLE WARUNKI PODPARCIA');
    disp(['WYZNACZNIK MACIERZ UK£ADU WYNOSI = ' num2str(d)]);
else
	x = A^-1 * B; % rozwi¹zanie uk³adu równañ ró¿nicowych
	w = x;
    % wyciêcie z wektora rozwi¹zañ ugiêæ we wêz³ach fikcyjnych
	w(n+4) = []; 
	w(n+3) = [];
	w(1) = [];
	w(1) = [];
    disp('ROZWI¥ZANIE - UGIÊCIE BELKI - WE WÊZ£ACH SIATKI');
	disp(w)
    wmax = max(w);
    wmin = min(w);
    figure;
    eps = 10^-10;
    wmax
    wmin
    if (abs(wmax-wmin)<eps)
        axis([0 L -(wmax+1) -(wmax-1)]);
    end
    p = plot([0 L],[0 0],'k-');
    set(p,'LineWidth',[3]);
    hold on
    grid on
	plot(0:h:L,-w,'ro-');
    title('UGIECIE BELKI');
    xlabel('x');
    ylabel('w(x)');
    legend('beam','w(x)');
end
