%
%                    PDE_code.m
%

  clear all

%
%  Set limits used to approximate forced heat equation
%

  L = 2;
  J = 40;
  K = 5;
  h = L/J;
  xj = 0:h:L;
  Nq = 11;
  hq = h/(Nq-1); 
  tf = 3;
  dt = 0.1;
  tvals = 0:dt:tf;
  
%
%  Set the parameter interval [pt1,pt2] and associted quadrature weights and points. Also
%  construct poins and weights for the standard interval [-1,1].  Finally construct the
%  normalization constants gamma.
%

  pt1 = 1.0;
  pt2 = 2.0;
  rho = 1/(pt2-pt1); 
  R = 15;
  hr = (pt2-pt1)/(R-1);
  wr = hr*ones(R,1);
  wr(1) = 0.5*hr;
  wr(R) = 0.5*hr;
  qr = pt1:hr:pt2;

%%
%  Input the Rs=10 Gauss-Legendre quadrature points and weights from a table.  Here qr_s 
%  denotes the points on the interval [-1,1] and qr_physical are points mapped to the 
%  physical interval [pt1,pt2] for the PDE evaluation.
%

  Rs = 10;
  qr_s = [-0.9739 -0.86506 -0.6794 -0.4334 -0.14887 0.14887 0.4334 0.6794 0.86506 0.9739];
  qr_physical = (1/2)*(pt1+pt2)*ones(size(qr_s)) + (1/2)*(pt2-pt1)*qr_s; 
  wr_s = 0.5*[0.06667 0.14945 0.21908 0.2693 0.2955 0.2955 0.2693 0.21908 0.14945 0.06667];

  gamma = [1/3 1/5 1/7 1/9 1/11];

%
% Compute MC samples for a direct evaluation.
%

  R_MC = 500;
  qr_MC = pt1 + (pt2-pt1)*rand(R_MC,1); 
  wr_MC = (1/R_MC)*ones(size(qr_MC));

%
%  Evaluate the Legendre polynomials at the quadrature points qr_s and construct products
%  with the weights.
%

  p1 = @(x) x;
  p2 = @(x) (3/2)*x.^2 - 1/2;
  p3 = @(x) (5/2)*x.^3 - (3/2)*x;
  p4 = @(x) (35/8)*x.^4 - (15/4)*x.^2 + 3/8;
  p5 = @(x) (63/8)*x.^5 - (70/8)*x.^3 + (15/8)*x; 
  P(1,:) = p1(qr_s);
  P(2,:) = p2(qr_s);
  P(3,:) = p3(qr_s);
  P(4,:) = p4(qr_s);
  P(5,:) = p5(qr_s);
  Pwt(1,:) = (1/gamma(1))*P(1,:).*wr_s;
  Pwt(2,:) = (1/gamma(2))*P(2,:).*wr_s;
  Pwt(3,:) = (1/gamma(3))*P(3,:).*wr_s;
  Pwt(4,:) = (1/gamma(4))*P(4,:).*wr_s;
  Pwt(5,:) = (1/gamma(5))*P(5,:).*wr_s;
  
%
% Create the right hand side vector
%

  xm1 = 0;
  xp1 = L;
%  gm = @(x) x.*((2-x).^4).*(x-xm1)/h;
%  gp = @(x) x.*((2-x).^4).*(xp1-x)/h;

  w = hq*ones(Nq,1); 
  w(1)  = 0.5*hq; 
  w(Nq) = 0.5*hq; 

  for j=1:J-1
    xm1 = xj(j+1)-h;
    xp1 = xj(j+1)+h;
    xm = xm1:hq:xj(j+1);
    xp = xj(j+1):hq:xp1;
    Gm = xm.*((xm-2).^4).*(xm-xm1)/h;
    Gp = xp.*((xp-2).^4).*(xp1-xp)/h;
    Im = Gm*w;
    Ip = Gp*w;
    g_vec(j,1) = Im+Ip;
  end

%
% Create the mass matrix M_mat, stiffness matrix K_mat, system matrix A_mat and vector G_vec
%

  M_mat = h*(diag((2/3)*ones(J-1,1)) + diag((1/6)*ones(J-2,1),-1) + diag((1/6)*ones(J-2,1),1));
  Minv = inv(M_mat);
  K_mat = (1/h)*(diag(2*ones(J-1,1)) + diag((-1)*ones(J-2,1),-1) + diag((-1)*ones(J-2,1),1));

%
% Create the vector of initial conditions
%

  r_init = sin(3*pi*xj(2:J)/2);

%
% Numerically integrate the ODE system for the collocation points qr
%

  rvals = zeros(R,J+1);
  for r = 1:R
    A_mat = -qr(r)*Minv*K_mat;
    G_vec = Minv*g_vec;
    ode_options = odeset('RelTol',1e-6);
    [t,Rvals] = ode45(@ODE_rhs,tvals,r_init,ode_options,A_mat,G_vec);
    rvals(r,2:J) = Rvals(length(t),:);
  end

%
% Numerically integrate the ODE system for the Legendre quadrature points qr_s
%

  rvals_s = zeros(Rs,J+1);
  for r = 1:Rs
    A_mat = -qr_physical(r)*Minv*K_mat;
    G_vec = Minv*g_vec;
    ode_options = odeset('RelTol',1e-6);
    [t,Rvals_s] = ode45(@ODE_rhs,tvals,r_init,ode_options,A_mat,G_vec);
    rvals_s(r,2:J) = Rvals_s(length(t),:);
  end

%
% Numerically integrate the ODE system using random points for MC
%

  rvals_MC = zeros(R_MC,J+1);
  for r = 1:R_MC
    A_mat = -qr_MC(r)*Minv*K_mat;
    G_vec = Minv*g_vec;
    ode_options = odeset('RelTol',1e-6);
    [t,Rvals_MC] = ode45(@ODE_rhs,tvals,r_init,ode_options,A_mat,G_vec);
    rvals_MC(r,2:J) = Rvals_MC(length(t),:);
  end

%
% Compute the mean and variance using collocation
%

  rmean_col = wr'*rvals*rho;
  var_col = zeros(1,J+1);
  for r=1:R
    var_col = var_col + wr(r)*rho*((rvals(r,:) - rmean_col).^2);
  end

  sigma_col = sqrt(var_col);

%
% Compute the mean and variance using discrete projection
%

  rmean_dp = wr_s*rvals_s;
  var_u_K = zeros(1,J+1); 
  for k=1:K
    u_K(k,:) = Pwt(k,:)*rvals_s;
    var_u_K = var_u_K + gamma(k)*u_K(k,:).^2; 
  end

%
% Compute the mean and variance using MC
%

  rmean_MC = wr_MC'*rvals_MC;
  var_MC = zeros(1,J+1);
  for r=1:R_MC
    var_MC = var_MC + wr_MC(r)*((rvals_MC(r,:) - rmean_MC).^2);
  end


%
% Plot the solution
%  

  figure(1)
  plot(xj,rmean_dp,'b-',xj,rmean_col,'r--',xj,rmean_col-2*sigma_col,'k:',xj,rmean_col+2*sigma_col,'k:','linewidth',3)
  axis([0 2 0 0.20]) 
  set(gca,'Fontsize',[20]);
  legend(' Discrete Projection', ' Collocation',' 2\sigma Interval','Location','South')
  xlabel('x-axis') 
  ylabel('Mean') 

  figure(2)
  plot(xj,var_u_K,'b-',xj,var_col,'--r',xj,var_MC,'k-.','linewidth',3)
  set(gca,'Fontsize',[20]);
  legend(' Discrete Projection',' Collocation',' Monte Carlo','Location','Northeast')
  xlabel('x-axis')
  ylabel('Variance')

  figure(3)
  contour(Rvals) 

  figure(4)
  plot(xj,var_MC,'k-.','linewidth',3)