%% Numerik D-MAVT, RW (Nipp), SS05
%
% This is the M-file solution to the ninth exercise sheet.
% To run, save this file to your account and ensure that the "Current directory"
% in MATLAB is the same as that where the file is saved.
%
%
% Then to run the file type "s9" in the command window


function s9
% Aufgabe 3
disp(['Aufgabe 3'])
x=[-0.90617985 -0.53846931 0          0.53846931 0.90617985]; % nodes 
w=[0.23692689   0.47862867 0.56888889 0.47862867 0.23692689]; % weights
I = pi/sqrt(2); % exact
% a)
Ia  = pi/2*sum(f(x,pi/2,pi/2).*w), errora=abs(I-Ia),
% b)
Ib1 = pi/6*sum(f(x,pi/6,pi/6).*w), 
Ib2 = pi/6*sum(f(x,pi/6,pi/2).*w),
Ib3 = pi/6*sum(f(x,pi/6,5*pi/6).*w),
Ib  = Ib1+Ib2+Ib3, errorb=abs(I-Ib),

% output:
% Ia =
%    2.23353745750014
% errora =
%    0.01209598842095
% Ib1 =
%    0.62655109238397
% Ib2 =
%    0.96834205198673
% Ib3 =
%    0.62655109238397
% Ib =
%    2.22144423675468
% errorb =
%    2.767675498560607e-06

%Aufgabe 4
disp(['Aufgabe 4'])
clear all
format long; format compact;
I=0.68299504214;	% exact solution
for TOL=[1e-1 1e-3 1e-6 1e-9]
  [Ia,rec] = adapquad(0,1,TOL,.5,'f2',0); % the algorithm
  Im = quad(@f2,0,1);		 % MATLAB quadrature
  disp(['TOL=',num2str(TOL),': |Ia-Im|=',num2str(abs(Ia-Im)),...
        ', |I-Ia|/I=',num2str(abs(I-Ia)/I),' Number of recursions=',num2str(rec)]);
end

%Aufgabe 4
%TOL=0.1: |Ia-Im|=0.0039365, |I-Ia|/I=0.0057633 Number of recursions=6
%TOL=0.001: |Ia-Im|=3.6489e-06, |I-Ia|/I=5.1163e-06 Number of recursions=27
%TOL=1e-06: |Ia-Im|=1.5423e-07, |I-Ia|/I=3.402e-10 Number of recursions=280
%TOL=1e-09: |Ia-Im|=1.5446e-07, |I-Ia|/I=1.6551e-12 Number of recursions=2916

% functions for aufgabe 3
function ff=f(x,a,b)
if nargin==1 a=1;b=0; end
ff=1./(1+cos(a*x+b).^2);
%---------------------------------------------------------
function df=dfexact(x,a,b)
if nargin==1 a=1;b=0; end
df=2*a*cos(a*x+b).*sin(a*x+b).*f(x,a,b).^2;
%---------------------------------------------------------

% functions for aufgabe 4
function [I,rec] = adapquad(a,b,er,Is,f,rec)

% adaptive quadrature to int_a^b f(x) dx \not = 0
% in:	a,b 	integrate over interval [a,b]
% 	er	realative tolerance
% 	Is    	estimate for the integral (\not = 0!)
%	f	the function to integrate (string)
% out:	I	approximation for int_a^b f(x) dx
%
%       Note Is need not be an accuarate for I, just should
%       be not equal to 0 and be of the correct order of
%       magnitude
%
%       Reference for algorithm is Schwartz pg 384-386     
% auth: fab@math.ethz.ch

h = b-a;                          
m = .5*(a+b);                
I1 = .5*h*(feval(f,a)+feval(f,b));	% trapezium rule
I2 = (I1+2*h*feval(f,m))/3.;            % Simpson rule
if (I1+Is*er/eps) == (I2+Is*er/eps)
 I = I2;
else % refine, recursion
 rec=rec+1;
 [Ia,rec] = adapquad(a,m,er,Is,f,rec);
 [Ib,rec] = adapquad(m,b,er,Is,f,rec);
 I=Ia+Ib;
end;
%---------------------------------------------------------------------

function ff = f2(x)
ff =  exp(-50*(x-.5).^2)+exp(-2*x);
