function MRSMES1D
close all
np = 33; % liczba wêz³ów
E = 2*10^8; % Modu³ Young'a
J = 0.1*(0.1^3)/12; % Pole przekroju
E = 1;
J = 1;
L = 2; % D³ugoœæ belki
obc = 1; % Schemat obci¹¿enia
ila = 40; % liczba siatek
wykres = 0;
nk = np + ila;

E1 = zeros(nk-np+1,4);
E2 = zeros(nk-np+1,4);
E3 = zeros(nk-np+1,4);
il = 0;
for n=np:2:nk
    h = L/(n-1);
    X = 0:h:L;
    A = zeros(n,n);
    B = zeros(n,1);
    for i = 2:n-1
        A(i,i-1:i+1) = [1 -2 1]/(h^2);
        B(i) = -mom(X(i),L,obc)/(E*J);
    end
    A(1,1) = 1;
    A(n,n) = 1;
    w1 = A^-1 * B;
    wp1 = zeros(n,1);
    for i=2:n-1
        wp1(i) = (w1(i+1) - w1(i-1))/(2*h);
    end
    wp1(1) = (w1(2) - w1(1))/h;
    wp1(n) = (w1(n) - w1(n-1))/h;
    w1sr = sqrt(dot(w1,w1)/n);
    w1max = max(abs(w1));
    wp1sr = sqrt(dot(wp1,wp1)/n);
    wp1max = max(abs(wp1));
    
    if nk==np
        subplot(2,1,1)
    end
    if nk>np
        subplot(2,2,1);
    end
    plot(X,w1,'ro-');
    hold on
    grid on
    view(0,-90);
    if nk==np
        subplot(2,1,2)
    end
    if nk>np
        subplot(2,2,3);
    end
    plot(X,wp1,'ro-');
    hold on
    grid on
    view(0,-90);
    
    N = n-1;
    A = zeros(n,n);
    B = zeros(n,1);
    for i=1:N
        K = [1 -1;-1 1]/h;
        f = zeros(2,1);
        Ng = 2;
        pg = [-1/sqrt(3) 1/sqrt(3)];
        for j=1:Ng
            PG = pg(j)*h/2 + h/2;
            f = f + Nk(PG,h) * (h/2)*mom(X(i)+PG,L,obc)/(E*J);
        end
        for j=1:2
            for k=1:2
                A(i-1+j,i-1+k) = A(i-1+j,i-1+k) + K(j,k);
            end
            B(i-1+j) = B(i-1+j) + f(j);
        end
    end
    A(n,:) = 0;
    A(1,:) = 0;
    A(:,n) = 0;
    A(:,1) = 0;
    B(1) = 0;
    B(n) = 0;
    A(1,1) = 1;
    A(n,n) = 1;
    w2 = A^-1 * B;
    wp2 = zeros(N,1);
    for i=1:N
        wp2(i) = (-1/h)*w2(i) + (1/h)*w2(i+1);
    end
    w2sr = sqrt(dot(w2,w2)/n);
    w2max = max(abs(w2));
    wp2sr = sqrt(dot(wp2,wp2)/N);
    wp2max = max(abs(wp2));
    
    if nk==np
        subplot(2,1,1)
    end
    if nk>np
        subplot(2,2,1);
    end
    plot(X,w2,'bo-');
    if nk==np
        subplot(2,1,2)
    end
    if nk>np
        subplot(2,2,3);
    end
    for i=1:N
        plot([X(i) X(i+1)],[wp2(i) wp2(i)],'bo-');
    end
    
    A = zeros(n,n);
    B = zeros(n,1);
    for i=1:N
        A(i,i) = A(i,i) + 1/h;
        A(i,i+1) = A(i,i+1) - 1/h;
        A(i+1,i) = A(i+1,i) - 1/h;
        A(i+1,i+1) = A(i+1,i+1) + 1/h;
        B(i) = B(i) + 0.5*h*mom(X(i),L,obc)/(E*J);
        B(i+1) = B(i+1) + 0.5*h*mom(X(i+1),L,obc)/(E*J);
    end
    A(n,:) = 0;
    A(1,:) = 0;
    A(:,n) = 0;
    A(:,1) = 0;
    B(1) = 0;
    B(n) = 0;
    A(1,1) = 1;
    A(n,n) = 1;
    w3 = A^-1 * B;
    wp3 = zeros(n,1);
    for i=2:n-1
        wp3(i) = (w3(i+1) - w3(i-1))/(2*h);
    end
    wp3(1) = (w3(2) - w3(1))/h;
    wp3(n) = (w3(n) - w3(n-1))/h;
    w3sr = sqrt(dot(w3,w3)/n);
    w3max = max(abs(w3));
    wp3sr = sqrt(dot(wp3,wp3)/n);
    wp3max = max(abs(wp3));
    
    
    if n>np
        il=il+1;
        e1sr = abs((w1sr-w1srpop)/w1sr);
        e2sr = abs((w2sr-w2srpop)/w2sr);
        e3sr = abs((w3sr-w3srpop)/w3sr);
        ep1sr = abs((wp1sr-wp1srpop)/wp1sr);
        ep2sr = abs((wp2sr-wp2srpop)/wp2sr);
        ep3sr = abs((wp3sr-wp3srpop)/wp3sr);
        e1max = abs((w1max-w1maxpop)/w1max);
        e2max = abs((w2max-w2maxpop)/w2max);
        e3max = abs((w3max-w3maxpop)/w3max);
        ep1max = abs((wp1max-wp1maxpop)/wp1max);
        ep2max = abs((wp2max-wp2maxpop)/wp2max);
        ep3max = abs((wp3max-wp3maxpop)/wp3max);
        E1(il,:) = [e1sr e1max ep1sr ep1max];
        E2(il,:) = [e2sr e2max ep2sr ep2max];
        E3(il,:) = [e3sr e3max ep3sr ep3max];
    end
    w1srpop = w1sr;
    w1maxpop = w1max;
    wp1srpop = wp1sr;
    wp1maxpop = wp1max;
    w2srpop = w2sr;
    w2maxpop = w2max;
    wp2srpop = wp2sr;
    wp2maxpop = wp2max;
    w3srpop = w3sr;
    w3maxpop = w3max;
    wp3srpop = wp3sr;
    wp3maxpop = wp3max;
    
    if nk==np
        subplot(2,1,1)
    end
    if nk>np
        subplot(2,2,1);
    end
    plot(X,w3,'go-');
    title(['UGIÊCIA: MRS LOK (W_m_a_x = ' num2str(max(w1),4) ') - MES (W_m_a_x = ' num2str(max(w2),4) ') - MRS WAR (W_m_a_x = ' num2str(max(w3),4) ')']);
    if nk==np
        subplot(2,1,2);
    end
    if nk>np
        subplot(2,2,3);
    end
    plot(X,wp3,'go-');
    title(['ODKSZTA£CENIA: MRS LOK (e_m_a_x = ' num2str(max(wp1),4) ') - MES (e_m_a_x = ' num2str(max(wp2),4) ') - MRS WAR (e_m_a_x = ' num2str(max(wp3),4) ')']);
    
    if nk>np & n>np
        kol = 'rbg';
        subplot(2,2,2);
        hold on
        grid on
        subplot(2,2,4);
        hold on
        grid on
        for k=1:3
            for i=2:2
                styl = '-';
                if i==2
                    styl = '--';
                end
                T = E1;
                if k==2
                    T = E2;
                end
                if k==3
                    T = E3;
                end
                for j=1:il
                    if abs(T(j,i))<10^-12
                        T(j,i) = 10^-12;
                    end
                    if abs(T(j,i+2))<10^-12
                        T(j,i+2) = 10^-12;
                    end
                end
                subplot(2,2,2);
                plot(log10(1:il),log10(T(1:il,i)),[kol(k) 'o' styl]);
                subplot(2,2,4);
                plot(log10(1:il),log10(T(1:il,i+2)),[kol(k) 'o' styl]);
            end
        end
    end
    if nk>np
        subplot(2,2,1);
        hold off
        subplot(2,2,3);
        hold off
        subplot(2,2,2);
        %l = legend('MRS^l^o^k_s_r','MRS^l^o^k_m_a_x','MES_s_r','MES_m_a_x','MRS^w^a^r_s_r','MRS^w^a^r_m_a_x');
        l = legend('MRS^l^o^k_m_a_x','MES_m_a_x','MRS^w^a^r_m_a_x');
        set(l,'FontSize',[8]);
        hold off
        subplot(2,2,4);
        hold off
    end
    if nk>np
        pause(0.001);
        if n==np
            pause;
        end
        if n==nk
            msgbox('KONIEC BADANIA ZBIEZNOSCI!','INFORMACJA','info','modal');
        end
    end
end
if (nk==np & wykres == 1)
    ile = 100;
    pause
    close all
    subplot(1,1,1);
    hold on
    grid on
    for i=1:N
        h1 = (X(i+1)-X(i))/(ile-1);
        r = zeros(ile,1);
        x = zeros(ile,1);
        pocz = i-2;
        kon = i+3;
        if kon>n
            kon = n;
        end
        if pocz<1
            pocz = 1;
        end
        il = kon - pocz + 1;
        A = zeros(il,4);
        B = zeros(il,1);
        W = zeros(il,il);
        for j=1:ile
            x(j) = X(i)+(j-1)*h1;
            for k=1:il
                d = X(pocz+k-1) - x(j);
                A(k,:) = [1 d 0.5*d^2 (1/6)*d^3];
                B(k) = w1(pocz+k-1);
                W(k,k) = 1/(d^6 + 10^-12);
            end
            a = (A'*W*A)^-1 * A'*W*B;
            r(j) = abs(a(3) + mom(x(j),L,obc)/(E*J));            
        end
        plot([X(i) x' X(i+1)],[0 r' 0],'r-');       
    end
    plot(X,zeros(n,1),'bo-');
end

function [z] = Nk(x,L)
z = [1-x/L;x/L];

function [z] = mom(x,L,obc)
if obc==1
    q = 1;
    z = 0.5*q*x*(L-x);
end
if obc==2
    P = 1;
    z = (P/2) * x;
    if x>L/2
        z = z - P*(x-L/2);
    end
end
if obc==3
    P = 1;
    z = (1/3)*P * x;
    if x>L/3
        z = z - P*(x-L/3);
    end
    if x>2*L/3
        z = z + P*(x-2*L/3);
    end
end
if obc==4
    q = 1;
    z = 0.5*q*x*(L-x);
    P = 1;
    z = z + (P/2) * x;
    if x>L/2
        z = z - P*(x-L/2);
    end
end
if obc==5
    z = sin(pi*x/L);
end
if obc==6
    z = sin(5*pi*x/L);
end