function func = numstab ( iexmpl );
% NumStab
%
% This m-file is part of the material of the lectures on
%   Numerical Mathematics
% at the Math Department of Halle University, Germany.
%
% Author :      Prof. Dr. M. Arnold, martin.arnold@mathematik.uni-halle.de
% Version of :  Dec 9, 2005
%
% Numerical (in)Stability in benchmark problems
%
% Parameters:
%   iexmpl (input)  : control flag "example"
%                       =1 ... (1-\cos x)/x
%                       =2 ... recursive evaluation of trigonometric
%                              functions
%
% Examples:
%   numstab ( 1 );        % print -dpng ../png/numstab1.png
%   numstab ( 2 );        % print -dpng ../png/numstab2.png

%
% -> Steuerparameter Ausgabe Bildschirm / PostScript
graphwidth  = 14.0;
graphheight = 14.0;
fac         =  1.0;
shift       =  0.5;
papershift  =  0.0;

ActFontSize   = 14;
ActLineWidth  =  2;
ActMarkerSize =  6;
RatioFlag     =  0;

border      = 1.0;

%
% -> Bestimme Groesse des Bildschirms und Skalierungsfaktor fuer die Ausgabe
root = 0;
set ( root, 'Units', 'centimeters' );
screensize = get ( root, 'ScreenSize' );
left   = screensize(1);
bottom = screensize(2);
width  = screensize(3);
height = screensize(4);
set ( root, 'Units', 'pixels' );
if RatioFlag==0,
  screenscale = 1.0;
elseif RatioFlag==1,
  screenscale = (width-2*border) / graphwidth;
end;

%
% -> Definiere den Bildschirmausschnitt
if length ( get ( root, 'Children' ) ) == 0,
  h = figure;
else
  h = gcf;
  clf;
end;
set ( h, 'Units',    'centimeters' );
set ( h, 'Position', ...
  [ left+width-screenscale*graphwidth-border height-2*border-screenscale*graphheight ...
    screenscale*graphwidth                   screenscale*graphheight ] );
set ( h, 'Units',    'pixels' );

%
% -> Skalierung der Ausgabe in PostScript
%
set ( gcf, 'PaperUnits',    'centimeters' );
set ( gcf, 'PaperPosition', [ papershift papershift graphwidth graphheight ] );

% -> plot
if iexmpl==1,

    % -> plot parameters
    logxmin = -12;
    logxmax = -4;
    fmin    = -6.0e-9;
    fmax    =  6.0e-9;

    % -> evaluate data
    nx  = 801;
    x   = logspace ( logxmin, logxmax, nx );
    f   = (1-cos(x))./x - x/2;
    res = x/2 .* ( -x.^2/12 + x.^4/360 - x.^6/(360*7*8) );

    % -> plot
    semilogx ( x, res, 'g', 'LineWidth', ActLineWidth );
    hold on;
    semilogx ( x, f, 'r' );
    hold off;
    xlabel ( 'x' );
    ylabel ( 'Ergebnis' );
    title ( '( 1 - cos x ) / x - x / 2' );
    axis ( [ 10^logxmin, 10^logxmax, fmin, fmax ] );

elseif iexmpl==2,

    % -> plot parameters
    mmin   = 0;
    mmax   = 15000;
    errmin = 1.0e-16;
    errmax = 1.0e-7;

    % -> evaluate data
    x = 1.0e-4;   
    cx = cos(x);   sx2q = sin(x/2)^2;   

    c3 = zeros ( mmax, 1 );
    c3(1) = 1.0;   c3(2) = cx;
    for i=3:mmax,
        c3(i) = 2*cx * c3(i-1) - c3(i-2);
    end;

    c2   = zeros ( mmax, 1 );
    cdel = zeros ( mmax, 1 );
    c2(1) = 1.0;   c2(2) = cx;   cdel(1) = -2*sx2q;
    for i=2:mmax-1,
        cdel(i) = -4*sx2q * c2(i) + cdel(i-1);
        c2(i+1) = c2(i) + cdel(i);
    end;

    % -> plot
    mvec  = (mmin:mmax-1)';
    cmvec = cos(mvec*x);
    semilogy ( ...
               mvec, abs((c3-cmvec)./cmvec), 'r'  ...
             , mvec, abs((c2-cmvec)./cmvec), 'g'  ...
             , 'LineWidth', ActLineWidth          ...
             );
    xlabel ( 'm' );
    ylabel ( 'Relativer Fehler' );
    title ( sprintf ( 'Rekursive Auswertung von  cos(mx)  fuer x = %6.1e', x ) );
    axis ( [ mmin mmax errmin errmax ] );

end;

%
% -> Korrektur der Groesse und Lage der Koordinatensysteme
allax = get ( gcf, 'Children' )';
for ax = allax,
  rect    = get ( ax, 'Position' );
  rect(2) = rect(2) + shift * (1-fac) * rect(4);
  rect(4) = fac * rect(4);
  set ( ax, 'Position', rect );
end;

%
% -> Aenderung der Schriftgroesse
for hf = allax,
  set ( hf, 'FontSize', [ActFontSize] );
  set ( get ( hf, 'title' ), 'FontSize', [ActFontSize] );
  set ( get ( hf, 'xlabel' ), 'FontSize', [ActFontSize] );
  set ( get ( hf, 'ylabel' ), 'FontSize', [ActFontSize] );
end;

% -> legend
if iexmpl==1,
    hl = legend ( 'Reihenentwicklung', 'numerisch' );
    for al = get ( hl, 'Children' )',
        for il=1:length(al),
            if ( strcmp(get(al(il),'Type'),'line')  ...
                 && sum(get(al(il),'Color')~=[1 0 0])==0 ),
                set ( al(il), 'LineWidth', ActLineWidth );
            end;
        end;
    end;
elseif iexmpl==2,
    legend ( 'Drei-Term-Rekursion'              ...
           , 'Gekoppelte Zwei-Term-Rekursionen' );
end;
