function ds = demo_spline_badcond()

    
%% Aproximovana funkce
%% a jeji druha derivace analyticky (jen pro vykresleni)
f      = inline('sin(x)','x');
fder2  = inline('-sin(x)','x');

%% Interval pro aproximaci
x_min = 0;
x_max = 2*pi;

%% Okrajove podminky: zadavame hodnoty derivace
yder_left  = 1.0;
yder_right = 1.0;
natural_left  = false;
natural_right = false;

%% Pocty bodu
bodu_input  = 50;   % Pocet "zadanych" bodu (vygeneruje se z f(x))
bodu_approx = 500;  % V kolika bodech chceme aproximovat
bodu_exact  = 500;  % Pocet bodu pro vykresleni analyt. funkce f''


%% Aproximace funkce z BODU_INPUT bodu na BODU_APPROX bodech pomoci splinu.

% vygenerujeme diskretni hodnoty
dx = (x_max-x_min) / (bodu_input-1);
x  = (x_min:dx:x_max);
y  = f(x);

% napocitame hodnoty druhych derivaci pro konstrukci splinu
% a) s pouzitim derivaci na obou okrajich
yder2 = NUMREC_spline( x , y , natural_left , yder_left, natural_right, yder_right );
% 2) s pouzitim prvni a druhe derivace vlevo 
%    (vlevo pouzijeme spoctenou y'' z predchoziho pripadu)
yder2_instab = zeros(size(yder2));
yder2_instab(1) = yder2(1);
yder2_instab(2) = 2*(  3* ( (y(2)-y(1))/dx - yder_left ) / dx ...
                     - yder2_instab(1) );
for i=3:bodu_input
    yder2_instab(i) =  6 * (y(i-2)-2*y(i-1)+y(i)) / dx^2 ...
                     - 4 * yder2_instab(i-1) - yder2_instab(i-2);
end;



% zvolime body a napocitame v nich presne a aproximovane hodnoty
dx_approx        = (x_max-x_min) / (bodu_approx-1);
x_approx         = (x_min:dx_approx:x_max);
y_exact          = f(x_approx);
y_approx         = zeros(size(x_approx));
y_approx_instab  = zeros(size(x_approx));
for i=1:bodu_approx
    y_approx(i)        = NUMREC_splint(x, y, yder2, x_approx(i));
    y_approx_instab(i) = NUMREC_splint(x, y, yder2_instab, x_approx(i));
end;

% stanovime rozmery obrazku dle velikosti obrazovky
scr = get( 0, 'Screensize' );
fig_wid = floor(scr(3)/3.5);
fig_off = floor(scr(3)/10);
fig_hei = floor(fig_wid/1.5);


% nakreslime
hf = figure;
set(hf,'Position',[fig_off fig_hei fig_wid fig_hei]);

ap = plot(x_approx,y_approx,'bo');
hold on;
api = plot(x_approx,y_approx_instab,'rx','linewidth',2);
hold off;
pl_max = max(y_approx(:));
pl_min = min(y_approx(:));
dpl = (pl_max-pl_min)/20;
axis([x_min x_max pl_min-dpl pl_max+dpl]);
title('f(x)','fontsize',16);
legend([ap api],{'y'' on both sides','y'' and y" at left'},'Location','northeast','fontsize',12)

% spocitame relativni chybu
err_y_approx        = abs((y_approx        - y_exact))./max(abs(y_exact),1e-14);
err_y_approx_instab = abs((y_approx_instab - y_exact))./max(abs(y_exact),1e-14);
err_y_approx        = log10(err_y_approx);
err_y_approx_instab = log10(err_y_approx_instab);

% nakreslime
hf = figure;
set(hf,'Position',[fig_off+fig_wid fig_hei fig_wid fig_hei]);
ap = plot(x_approx,err_y_approx,'bo');
hold on;
api = plot(x_approx,err_y_approx_instab,'rx','linewidth',2);
hold off;
pl_max = max(err_y_approx_instab(:));
pl_min = min(err_y_approx_instab(:));
dpl = (pl_max-pl_min)/20;
axis([x_min x_max pl_min-dpl pl_max+dpl]);
title('{log_{10}} [ {L_{1}^{rel}}{(y_i)} ]','fontsize',16);
legend([ap api],{'y'' on both sides','y'' and y" at left'},'Location','northwest','fontsize',12)



%% Vykresleni druhych derivaci pouzitych ke konstrukci splinu
%% vypoctenych obema zpusoby

% napocitame hustou sit s presnou (analytickou) druhou derivaci
dx_exact    = (x_max-x_min) / (bodu_exact-1);
x_exact     = (x_min:dx_exact:x_max);
yder2_exact = fder2(x_exact);

% nakreslime je spolu s hodnotami druhe derivace spoctenymi pri
% konstrukci splinu
hf = figure;
set(hf,'Position',[fig_off+2*fig_wid fig_hei fig_wid fig_hei]);

ex = plot(x_exact,yder2_exact,'k-');
hold on;
ap  = plot(x,yder2,'bo');
api = plot(x,yder2_instab,'rx','linewidth',2);
hold off;
pl_max = max(yder2_exact(:));
pl_min = min(yder2_exact(:));
dpl = (pl_max-pl_min)/20;
axis([x_min x_max pl_min-dpl pl_max+dpl]);
title('f"(x)','fontsize',16);
legend([ex, ap, api],{'exact','y'' on both sides','y'' and y" at left'},...
       'Location','southeast','fontsize',12)



end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Subrutiny pro konstrukci a vypocet splinu %
%         podle Numerical Recipes           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%% Subrutina pro vypocet druhych derivaci
function [yder2] = NUMREC_spline(x, y, l_natur, l_yder, r_natur, r_yder)
   % Calculate second derivatives for the interpolation by cubic spline
   %
   % Written according to routine SPLINE from Numerical Recipes
   %
   % INPUT:
   % - Arrays X and Y of length NPT containing a tabulated function, 
   %   i.e., y_i = f(xi), with x_1 < x_2 < .. . < x_NPT
   % - On each boundary either flag to use natural spline
   %   or value of first derivative      
   %    
   % OUTPUT: 
   % - array YDER2 of length npt which contains the second derivatives
   %   of the interpolating function at the tabulated points X_i. 
       
   npt   = length(x);
   u     = zeros(size(x));
   yder2 = zeros(size(x));
   %
   if (l_natur)  % Left BC natural
       yder2(1) = 0.0;
       u(1)     = 0.0;
   else          % Left BC specified by first derivative l_yder
       yder2(1)  = -0.5;
       u(1)      = (3.0/(x(2)-x(1))) * ((y(2)-y(1)) / (x(2)-x(1)) - l_yder);
   end;
   % Decomposition loop of the tridiagonal algorithm. 
   % yder2 and u are used for temporary storage of the decomposed factors.
   for i = 2:npt-1
       sig      = (x(i)-x(i-1)) / (x(i+1)-x(i-1));
       p        = sig * yder2(i-1) + 2.0;
       yder2(i) = (sig-1.0) / p;
       u(i)     = ( 6. * (  (y(i+1)-y(i)) / (x(i+1)-x(i)) ...
                           -(y(i)-y(i-1)) / (x(i)-x(i-1)) ) ...
                       / (x(i+1)-x(i-1)) - sig*u(i-1) ...
                  ) / p;
   end;
   %
   if (r_natur)  % Right BC natural
       qn = 0.0;
       un = 0.0;
   else          % Right BC specified by first derivative r_yder
       qn = 0.5;
       un = (3.0/(x(npt)-x(npt-1))) ...
            * (r_yder - (y(npt)-y(npt-1))/(x(npt)-x(npt-1)));
   end;
   % Backsubstitution loop of the tridiagonal algorithm.
   yder2(npt) = (un-qn*u(npt-1)) / (qn*yder2(npt-1) + 1.0);
   for k=npt-1:-1:1
      yder2(k) = yder2(k) * yder2(k+1) + u(k);
   end;
end


%%% Subrutina pro vypocet interpolovane hodnoty
function y = NUMREC_splint(xa, ya, yader2, x)
   % Calculate the interpolated value at point x by cubic spline
   %
   % Written according to routine SPLINT from Numerical Recipes
   %
   % INPUT:
   % - Arrays XA and YA of length NPT containing a tabulated function, 
   %   i.e., y_i = f(xi), with x_1 < x_2 < .. . < x_NPT
   % - Arrays XADER2 containing the second derivatives
   % - Point X
   %    
   % OUTPUT: 
   % - Value of Y(X) approximated by cubic spline

   npt   = length(xa);
   %    
   % We will find the right place in the table by means of bisection.
   % This is optimal if sequential calls to this routine are at random
   % values of x. If sequential calls are in order, and closely spaced,
   % one would do better to store previous values of klo and khi
   % and test if they remain appropriate on the next call.
   klo = 1;
   khi = npt;
   while (khi-klo > 1)
       k = floor ((khi+klo)/2);
       if (xa(k) > x) 
           khi = k; 
       else
           klo = k;
       end;
   end;
   % klo and khi now bracket the input value of x.
   h = xa(khi)-xa(klo);
   if (h==0.0)
       disp(['Error in routine NUMREC_splint. Points of XA must be distinct.']);
       return;
   end; 
   % Cubic spline polynomial is now evaluated.
   a = (xa(khi)-x) / h; 
   b = (x-xa(klo)) / h;
   y = a*ya(klo) + b*ya(khi) ...
       + ( (a^3-a)*yader2(klo) + (b^3-b)*yader2(khi) ) * h^2 / 6.0;
end


