function u = reseni_ukol1(k,l,c,n)
% u = reseni_ukol1(k,l,c,n)
%
% Medota siti pro danou 1D Dirichletovu ulohu s konstantni 
% hustotou zatezujici sily
%
% k   ... tuhost struny
% l   ... delka struny
% c   ... konstantni hustota zatezujici sily 
% n+1 ... pocet pravidelnych uzlu na strune
%
% u   ... vysledny pruhyb struny 


% vzdalenost mezi sousednimi uzly
h = l / n;


% vytvoreni matice tuhosti K
K = diag(2*ones(n-1,1)) + diag(-ones(n-2,1),1) + diag(-ones(n-2,1),-1);
K = k/h^2 * K;


% vektor prave strany f
f = c*ones(n-1,1);


% reseni soustavy u
u = K \ f;
u = [0; u; 0];


% vykresleni reseni
figure; hold on
x = 0 : h : l;
plot(x,u)


% srovnani s analytickym resenim
u_a = c/2/k*x.*(l-x);
plot(x,u_a,'r')
% maximalni absolutni chyba
err_abs = max(abs(u - u_a')); 
% relativni chyba
err_rel = norm(u - u_a') / norm(u_a);
% vypis vysledku
disp(['Maximalni absolutni chyba: ', num2str(err_abs)])
disp(['Relativni chyba: ', num2str(err_rel)])

% konec

