clear all N = 5; beta = 10; mu = 1.0; sigma = 0.25; sigma_beta = 2; time = 5; dt = .1; Ntime = time/dt; % % Define the time grid and Hermite quadrature points and weights. The vectors alpha_vals % and beta_vals contain the associated parameter values and u contains the state values. % tvals = [0:dt:time]; zvals = sqrt(2)*[-1.6507 -.5246 .5246 1.6507]; weights = [8.1313e-2 8.0491e-1 8.0491e-1 8.1313e-2]; alpha_vals = mu*ones(size(zvals)) + zvals*sigma; beta_vals = beta*ones(size(zvals)) + zvals*sigma_beta; M = length(zvals); for j1 = 1:M for j2 = 1:M u(j1,j2,:) = beta_vals(j2)*exp(-alpha_vals(j1)*tvals); end end % % Construct the normalization constants and Hermite polynomials. Initialize % the cofficient vectors. % h0 = factorial(0); poly0 = ones(1,M); h1 = factorial(1); poly1 = zvals; h2 = factorial(2); poly2 = zvals.^2 - ones(1,M); w0 = zeros(1,Ntime+1); w1 = zeros(1,Ntime+1); w2 = zeros(1,Ntime+1); w3 = zeros(1,Ntime+1); w4 = zeros(1,Ntime+1); w5 = zeros(1,Ntime+1); % % Use Gauss-Hermite quadrature to construct the coefficient vectors for each % time value. % for k = 1:Ntime+1 for j1 = 1:M for j2 = 1:M w0(k) = w0(k) + (1/(pi*h0*h0))*u(j1,j2,k)*poly0(j1)*poly0(j2)*weights(j1)*weights(j2); w1(k) = w1(k) + (1/(pi*h1*h0))*u(j1,j2,k)*poly1(j1)*poly0(j2)*weights(j1)*weights(j2); w2(k) = w2(k) + (1/(pi*h0*h1))*u(j1,j2,k)*poly0(j1)*poly1(j2)*weights(j1)*weights(j2); w3(k) = w3(k) + (1/(pi*h2*h0))*u(j1,j2,k)*poly2(j1)*poly0(j2)*weights(j1)*weights(j2); w4(k) = w4(k) + (1/(pi*h1*h1))*u(j1,j2,k)*poly1(j1)*poly1(j2)*weights(j1)*weights(j2); w5(k) = w5(k) + (1/(pi*h0*h2))*u(j1,j2,k)*poly0(j1)*poly2(j2)*weights(j1)*weights(j2); end end end W = [w0; w1; w2; w3; w4; w5]; gamma = [h0*h0 h1*h0 h0*h1 h2*h0 h1*h1 h0*h2]; % % Construct the variance values. % for n = 1:Ntime+1 var_sum = 0; for j = 2:N+1 var_sum = var_sum + (W(j,n)^2)*gamma(j); end var(n) = var_sum; end % % Construct the approximate and true mean, variance and 2 sigma % confidence values. % mean = W(1,:); sd = sqrt(var); mean2p = mean + 2*sd; mean2m = mean - 2*sd; det = exp(-mu*tvals)*beta; mean_true = beta*exp(-mu*tvals).*exp(.5*(sigma^2)*(tvals.^2)); var_true = exp(-2*mu*tvals).*(exp(2*(sigma^2)*(tvals.^2))*(beta^2 + sigma_beta^2) - exp((sigma^2)*(tvals.^2))*beta^2); sd_true = sqrt(var_true); mean2p_true = mean_true + 2*sd_true; mean2m_true = mean_true - 2*sd_true; figure(1) plot(tvals,mean,tvals,det,'r',tvals,mean_true,'-g',tvals,mean2p,'--k',tvals,mean2m,'--k','linewidth',3) hold on plot(tvals,mean2p_true,'--g',tvals,mean2m_true,'--g','linewidth',3) hold off set(gca,'Fontsize',[20]); xlabel('Time (s)') ylabel('Displacement (m)') legend('Mean','Deterministic','True Mean','Location','NorthEast')