
a=0;
b=2;
syms m
%m=(a+b)/2;
r=(b-a)/2;
uI=1;

% Jacobi matrix
N=2;
ia=1:N;
ib=0:N;
ic=0:N-1;
Ma=ia./(2*ia+1);
Mb=0*ib;
Mc=(ic+1)./(2*ic+1);
M=diag(Ma,1)+diag(Mb)+diag(Mc,-1);
disp('Stochastic Galerkin results in modal system du/dt=Au:')
A=subs(r*M+m*eye(N+1),m,1)

% Eigenvalue decomposition of A
[V,D]=eig(A);
disp(['Eigenvalues of A are '])
disp(diag(D)')
disp(['which expand to ', num2str(double(diag(D))')]) % show decimal

disp(sprintf('\nFigure shows the numerical solution of modal equations'))
disp('where var(u)=u_1^2+u_2^2')
A=double(A); % make A not symbolic for numerical simulations
[t,u]=ode45(@(t,u) -A*u,[0 10],[uI,0,0]);
figure
plot(t,u(:,1),'-+r');
hold on
plot(t,u(:,1)+2*sqrt(u(:,2).^2+u(:,3).^2),'-.k')
plot(t,u(:,1)-2*sqrt(u(:,2).^2+u(:,3).^2),'--b')
xlabel('x')
ylabel('Solution to y''=-kx, k ~ U[0,2]')
%print('-dpng','examplePSCS')


disp(sprintf('\nEigenvalue approach to solving modal equations involves change'))
disp('of variables including')
wI=inv(V)*[uI,0,0]'
disp('and the matrix exponential given by')
syms s
expDt=diag(exp(-diag(D*s)))
disp('converting back gives, which matches the plot from ode45')
uN=simplify(V*expDt*wI)
disp('note that the weights for u0 are thus')
double(V(1,:)'.*wI)'
plot(t,subs(uN(1),s,t),'-xg') % should match ode45 solution
legend('u_0 (SG via ode45)','u_0+2\sigma','u_0-2\sigma','u_0 (SG analytical)')

disp(sprintf('\nStochastic collocation involves sampling the ode at collocation nodes'))
% Roots of 3rd degree Legendre poly gives the following collocation nodes
syms x
roots = vpasolve(legendreP(3,x) == 0)';
nodes=subs(r*roots+m,m,1);
disp(num2str(double(nodes)));

disp('and taking a weighted average using weights')
wi=[5/9,8/9,5/9]/2 
%note that the extra factor 1/2 is due to the fact that dF = 1/2 dx
%standard Guassian quadrature tables assume weight w(x)=1.
disp('which is identical to the eigenvalue decomposition of SG')
