//
// UQ avec méthode du chaos polynomial, toy problem
//
///  3 direct computations
//
amin=1;
amax=3;
  a=-sqrt(3/5)+2
  N=10000
  S1=1
  for i=1:N
    S1=S1+1/N*(-a*S1) 
  end
///////////////////////////
  a=0+2
  N=10000
  S2=1
  for i=1:N
    S2=S2+1/N*(-a*S2) 
  end
///////////////////////////
  a=sqrt(3/5)+2
  N=10000
  S3=1
  for i=1:N
    S3=S3+1/N*(-a*S3) 
  end
////////////////////////////////////
deff('y=L0(x)','y=1')
deff('y=L1(x)','y=x')
deff('y=L2(x)','y=3/2*x^2-1/2')
//
alpha0=5/9*S1*L0(-sqrt(3/5))+8/9*S2*L0(0)+5/9*S3*L0(sqrt(3/5));
alpha1=5/9*S1*L1(-sqrt(3/5))+8/9*S2*L1(0)+5/9*S3*L1(sqrt(3/5));
alpha2=5/9*S1*L2(-sqrt(3/5))+8/9*S2*L2(0)+5/9*S3*L2(sqrt(3/5));
//
//
disp('exact mean of S');disp(1/2*(exp(-1)-exp(-3)));
disp('mean of S with polynomial chaos');disp(alpha0/2);
//
disp('exact variance of S');disp(1/4*(%e^(-2)-%e^(-6))-(1/2*(exp(-1)-exp(-3)))^2)
varappro=1/2*((alpha1*1/sqrt(2/3))^2+(alpha2*1/sqrt(2/5))^2)
disp('variance of S with polynomial chaos');disp(varappro)
