%lab7.m computing integrals in octave....

%a few simple integrals using different quad fucntions:
%quad(), quadv() quadl() quadgk() quadcc()
%let compute \int_0^1 x^p
int=e=zeros(5,1);
format short e
for p=2:4:20,
int(1)=quad(f=@(x) x.^p,0,1);
e(1)=int(1)-1/(p+1);

int(2)=quadv(f,0,1);
e(2)=int(2)-1/(p+1);

int(3)=quadl(f,0,1);
e(3)=int(3)-1/(p+1);

int(4)=quadgk(f,0,1);
e(4)=int(4)-1/(p+1);

int(5)=quadcc(f,0,1);
e(5)=int(5)-1/(p+1);
[p*ones(size(int)),int,e]
end

%singular functions (not all above quad?? functions works then...)
%e.g. quadl() does not (why?)

  p=-0.5;
int(1)=quad(f=@(x) x.^p,0,1);
e(1)=int(1)-1/(p+1);

int(2)=quadv(f,0,1);
e(2)=int(2)-1/(p+1);

int(3)=-1;%quadl(f,0,1);
e(3)=-1;%int(3)-1/(p+1);

int(4)=quadgk(f,0,1);
e(4)=int(4)-1/(p+1);

int(5)=quadcc(f,0,1);
e(5)=int(5)-1/(p+1);
[p*ones(size(int)),int,e]

%we can try to use quadl() by giving a different integrating interval i.e
%this integral approx. the original one

int(3)=quadl(f,1e-20,1)
e(3)=int(3)-1/(p+1)
%quadl() is not he best choice for a integrand with a singularity at an end...


%let's try infinite interval just with the quad() one
int=quad(f=@(x) exp(-x^2/2),0,inf)
int-sqrt(0.5*pi)


%plot the graph of F(s)=\int_\infty^s exp(-x^2/2) dx on [-20,20]
%simple straightforward solution
f=@(x) exp(-x.*x*0.5); %integrand
F=@(s) quad(f,-inf,s);
s=linspace(-20,20);
tic;
for k=1:length(s),
  Fs(k)=F(s(k));
end
t1=toc
plot(s,Fs);pause(2);
%we want to solve (F(s^*)=0.5
ss=fzero(@(s) F(s)-0.5,[-10,10])
F(ss)
%not that F(s+h)=F(s)+\int_s^(s+h) f(z)dz
%and this integral should be easier (and sometimes) cheaper to compute 
%naturally unless a loop usage which would slow down the process
%however if a standard call to quad() fails one can try this simple trick
tic;
Fs1(1)=F(s(1));
quad_options("absolute tolerance",1.5e-9);
quad_options("relative tolerance",1e-8)
for k=2:length(s),
  Fs1(k)=Fs1(k-1)+ quad(f,s(k-1),s(k));
end
t2=toc

plot(s,Fs-Fs1);pause(2);

%L^2 norm
function y=Ltwonorm(f,a,b)
  y=sqrt(quad(@(x) f(x)*f(x),a,b));
  
end  
%let's test for f(x)=|x| on [-1,1]
l2n=Ltwonorm(@(x) abs(x),-1,1)
l2n-sqrt(2/3)


%L^p norm
function y=Lpnorm(f,a,b,p=2)
  if(p>=1)
   y=(quad(@(x) abs(f(x))^p,a,b))^(1/p);
  else
    warning("p should be >=1.");
    return;
  end
end  
%let's test for f(x)=|x| on [-1,1]
l2n=Lpnorm(@(x) abs(x),-1,1)
l2n-sqrt(2/3)

p=3;
l3n=Lpnorm(@(x) abs(x),-1,1,p)
l3n-(2/(p+1))^(1/p)
%try different p...

%try to approximate W^{1,2} norm i.e.
%||f||=sqrt(||f||_{L^2} + ||f'||_{L^2})

%best approximation problem in  L^2[a,b]
%find a polynomial w of degree <=p s.t. 
%||f-w||_{L^2(a,b)} minimal...
% let assume that we want to find coeff in power basis
% the nwe have to solve the system Hx=F with
%H athe Hilbr matrix (hilb(p+1)) and F=(\int_a^b fx^k dx)_{k=0,..,p}

function [w,er]=L2bapoly(f,a,b,p)
%we compute bet approximation poly of degree <=p
%input 
%f function handle to math function : y=f(x)
%a,b - the inteval endswitch
%p - degree of polynomial space
%OUTPUT
%w - a vector of  poly coeff (in octave standard i.e. in the basis (x^p,...,x^0)
%er - L^2 error ||f- \sum_(k=1)^(p+1) w(k)x^(p+1-k)||_{L^2(a,b)}
f1=@(t) f(a+(b-a)*t);
H=hilb(p+1);
F=zeros(p+1,1);
for k=0:p,
  F(k+1)=quad(@(x) x.^k.*f1(x),0,1);
end  
F
%F=flipud(F);
w=H\F;
w=flipud(w);
ab201=@(x) (x-a)/(b-a);
g=@(x) (f(x)-polyval(w,ab201(x)));
er=quad(@(x) g(x).*g(x),a,b);
er=sqrt(er);
endfunction

f=@sin;
a=0;b=1;
p=3;
[w,er]=L2bapoly(f,a,b,p)
x=linspace(a,b);
plot(x,fx=f(x),x,px=polyval(w,(x-a)/b-a));pause(2);
plot(x,fx-px);pause(2);

a=-1;b=4;
f=@(x) x.*x-2;
[w,er]=L2bapoly(f,a,b,p)
x=linspace(a,b);
plot(x,fx=f(x),x,px=polyval(w,(x-a)/(b-a)));pause(2);
plot(x,fx-px);pause(2);



a=-1;b=1.5;p=8;
f=@(x) (x>0); %characteristic function of (0,+\infty)
[w,er]=L2bapoly(f,a,b,p)
x=linspace(a,b);
plot(x,fx=f(x),x,px=polyval(w,(x-a)/(b-a)));pause(2);
plot(x,fx-px);pause(2);



a=-1;b=1.5;p=8;
f=@(x) abs(x); %|x|
[w,er]=L2bapoly(f,a,b,p)
x=linspace(a,b);
plot(x,fx=f(x),x,px=polyval(w,(x-a)/(b-a)));pause(2);
plot(x,fx-px);pause(2);


%double integral
%integrate \int_({x^2+y^2<=1} 1 dx dy
%it is enough to do it for x,y>0
yx=@(x) sqrt(1-x.*x);
F=@(x,y) 1.0*(x.*x+y*y-1<0)

i1=dblquad(F,0,1,0,1)
%using quadv() twice:
%\int_0^1 \int_0^(1-x^2) dy  dx
i2=quadv(@(x) quadv(@(x) ones(size(x)),0,sqrt(1-x*x)),0,1)
i2-i1
i=pi/4
0.25*pi-[i1,i2]
%so 1st method better...try a different quad??() for the 2nd one...

%trapz() function - if we have certain valus of f at given x
%(obtained by some computations/or experiments/measuremts etc
%we can approximated an integral using trapz() i.e. composite trapezoid rule
N=51;
x=linspace(0,1,N);
x2=x.*x;
int=trapz(x,x2)
er(1)=int-1/3 %1e-5 orderof error
%take more points
x=linspace(0,1,2*(N-1));
x2=x.*x;
int=trapz(x,x.*x)
er(2)=int-1/3 
er(1)/er(2)
%we can check experimentally the order 
N=20;M=7;
for k=1:M,
  x=linspace(0,1,N+1);
  int=trapz(x,x.*x);
  er(k)=int-1/3;
  N*=2;
end
er
er(1:M-1)./er(2:M)
%results suggests that error behaves like O(N^(-2))
%please repeat it for e.g. sin(x) or x^(k+0.5) k=0,1,2