We are considering the spring model:
\begin{align*}
\frac{d^2z}{dt^2} + kz &= 0 \\[10pt]
z(0) = 3, \quad \frac{dz}{dt}(0) &= 0
\end{align*}
\)
Which has a solution given by:
\begin{align*}
z(t,k) = 3\cos(\sqrt{k}\cdot t)
\end{align*}
\)
We construct a Legendre surrogate \(f_S^K(t,\xi)=\sum\limits_{k=0}^{K}u_k\Psi_k(t)\) using the \(K+1=4\) basis functions:
\begin{align*}
\Psi_0(q)=1, \quad \Psi_1(q)=q, \quad \Psi_2(q)=\frac{1}{2}(3q^2-1), \quad \Psi_3(q)=\frac{1}{2}(5q^3-3q)
\end{align*}
\)
These are the first four Legendre polynomials. They are orthogonal on \(\Gamma=[-1,1]\) with respect to the uniform distribution. The weights \(u_k\) are found via discrete projection; \(u_k=\frac{1}{\gamma_k}\int_{\Gamma}f(t)\Psi_k(t)\rho(t) \, dt=\frac{1}{2k+1}\int_{\left[-1,1\right]}\frac{1}{2}z(t)\Psi_k(t) \, dt\). We take \(k \sim \mathcal{U}(\bar{k} – \sigma_k, \bar{k} + \sigma_k)\) and \(\xi \in \left[-1,1\right]\) with \(\bar{k} = 8.5\) and \(\sigma_k = 0.001\).
In addition to the approximation with \(R=10\) Legendre quadrature points, we use \(M=10^5\) samples to perform Monte Carlo approximations. The comparison of the means and standard deviations are shown in the figures below.
%%%2a. set up%%%
K=3; %Number of tensored Legendre basis functions Psi(q)%
N_MC=1e5; %Number of Monte Carlo samples%
M=10; %Number of Gauss-Legendre quadrature points used to approximate discrete projection%
mu_k=8.5; %k \in U(kbar-sigma_k, kbar+sigmak)%
sigma_k=0.001; %k \in U(kbar-sigma_k, kbar+sigmak)%
%%%2b. Compute Monte Carlo Samples%%%
q=(mu_k-sigma_k)+2*sigma_k*rand(N_MC,1);
pt1=0; %lower time%
pt2=5; %upper time%
dt=0.005; %change in time%
t=[pt1:dt:pt2];
N_t=length(t);
%%%2c. Input the M=10 Gauss-Legendre quadrature points and weights from a table%%%
%Recall that one must multiply the weights by 0.5 to account for the density \rho(q) = 1/2%
zvals=[-0.9739 -0.86506 -0.6794 -0.4334 -0.14887 0.14887 0.4334 0.6794 0.86506 0.9739];
weights=0.5*[0.06667 0.14945 0.21908 0.2693 0.2955 0.2955 0.2693 0.21908 0.14945 0.06667];
sigmak = 2*sigma_k/sqrt(12);
kvals = mu_k*ones(size(zvals)) + sqrt(3)*sigmak*zvals;
%%%2d. Constuct the Legendre polynomials P_0, P_1, and P_2 evaluated at the quadrature points zvals%%%
%Also construct the coefficients h0, h1, and h2 used to construct gamma_k%
poly0=ones(1,M); %P_0(q)=1%
poly1=zvals; %P_1(q)=q%
poly2=(3/2)*zvals.^2-(1/2)*ones(1,M); %P_2(q)=(1/2)(3q^2-1)%
poly3=(5/2)*zvals.^3-(3/2)*zvals; %P_3(q)=(1/2)(5q^3-3q)%
h0=1; %gamma_k=1/(2k+1)%
h1=1/3;
h2=1/5;
h3=1/7;
y_k0=zeros(1,N_t); %Initialize the K+1 Legendre polynomials%
y_k1=zeros(1,N_t); %Initialize the K+1 Legendre polynomials%
y_k2=zeros(1,N_t); %Initialize the K+1 Legendre polynomials%
y_k3=zeros(1,N_t); %Initialize the K+1 Legendre polynomials%
%%%2e. Compute coefficients%%%
for k=1:N_t
deterministic(k)=3*cos(sqrt(mu_k)*t(k));
MC=3*cos(sqrt(q)*t(k));
MC_mean(k)=(1/N_MC)*sum(MC);
tmp=0;
for j=1:N_MC
tmp=tmp+(MC(j)-MC_mean(k))^2;
end
MC_sigma(k)=sqrt((1/(N_MC-1))*tmp);
for j=1:M
yval=3*cos(sqrt(kvals(j))*t(k)); %high fidelity%
y_k0(k)=y_k0(k)+(1/h0)*yval*poly0(j)*weights(j);
y_k1(k)=y_k1(k)+(1/h1)*yval*poly1(j)*weights(j);
y_k2(k)=y_k2(k)+(1/h2)*yval*poly2(j)*weights(j);
y_k3(k)=y_k3(k)+(1/h3)*yval*poly3(j)*weights(j);
end
end
%%%2f. Compute mean and variance%%%
Y_k=[y_k0; y_k1; y_k2; y_k3];
gamma=[h0 h1 h2 h3];
for n=1:N_t
var_sum=0;
for j=2:K+1
var_sum=var_sum+(Y_k(j,n)^2)*gamma(j);
end
var(n)=var_sum;
end
DP_mean=Y_k(1,:); %Discrete projection mean%
DP_sd=sqrt(var); %Discrete projection sd%
%%%2g. Plot%%%
figure(3)
plot(t,DP_mean,'k-',t,MC_mean,'--b','linewidth',3)
set(gca,'Fontsize',[22]);
xlabel('t')
ylabel('Response Mean')
legend('Discrete Projection','Monte Carlo','Location','NorthEast')
figure(4)
plot(t,DP_sd,'-k',t,MC_sigma,'--b','linewidth',3)
set(gca,'Fontsize',[22]);
xlabel('t')
ylabel('Response Standard Deviation')
legend('Discrete Projection','Monte Carlo','Location','NorthWest')