% Source code for the figure for final exam question 1
%
% There is a small cheat: To avoid that all the graphs lie on top
% of each other, some of the functions are scaled, resulting in a
% vertical shift of the graph on the log-log error plot.

N = 2.^[2:16];

function S = trapez (fun,a,b,N)
  h = (b-a)/N;
  fy = feval(fun,linspace(a,b,N+1));
  fy(1) = fy(1)/2;
  fy(N+1) = fy(N+1)/2;
  S = h*sum(fy);
end

function y = f(x)
  y = (sin(x)).^2;
end

function y = g(x)
  y = 1/100*x.^2;
end

function y = h(x)
  y = 10*sqrt(x);
end

function y=j(x)
  y = (x<1);
end

for k=1:length(N);
  f_err(k) = abs(pi-trapez('f',0,2*pi,N(k)));
  g_err(k) = abs(1/100*8/3*pi^3-trapez('g',0,2*pi,N(k)));
  h_err(k) = abs(10*4/3*sqrt(2)*pi^(3/2)-trapez('h',0,2*pi,N(k)));
  j_err(k) = abs(1-trapez('j',0,2*pi,N(k)));
end

title('Trapezoidal Rule');
xlabel('N');
ylabel('Error');
loglog(N, 1./(N.^2), ';Slope 2;', ...
       N, 30./(N.^1.5), ';Slope 3/2;', ...
       N, f_err, 'x;f(x);', ...
       N, g_err, 'o;g(x);', ...
       N, h_err, '*;h(x);', ...
       N, j_err, '+;j(x);')

gset term postscript eps monochrome;
gset output "/tmp/loglog.eps"
replot
gset term x11