# this works for only one random variable 
# and the polynomials must be orthonormal
param N; # size of matrix and vector
param M; # number of elements considered in the Fourier Expansion of the distribution
param L{1..M,1..M}; # <P_i, zeta P_j>
param LL{1..M,1..M,1..M}; # <P_i, P_j  P_k> 
param Q {i in 1..N, j in 1..N} 
	= if (i+j < 3) then 1
                                          else if (i+j > (2*N-1)) then 1
		else if (i=j) then 2
		                    else if (abs(i-j)=1) then -1;

param DQ {i in 1..N, j in 1..N} 
	= if (i = j) then 2*(i/N)*(1-i/N)*cos(i/N) else 0; 

param RHS{1..M}; # <P_i>
# total matrix is Q+Uniform(-1,1)*DQ

var x{1..M,1..N}; 

minimize Obj: sum{i1 in 1..M, j1 in 1..N, j2 in 1..N} (x[i1,j1]*Q[j1,j2]*x[i1,j2])+sum{i1 in 1..M, i2 in 1..M, j1 in 1..N, j2 in 1..N} (x[i1,j1]*L[i1,i2]*DQ[j1,j2]*x[i2,j2]);

# Galerkin Constraints 

subject to galerkin {i1 in 1..M}: sum{i2 in 1..M, i3 in 1..M, j2 in 1..N} (x[i2,j2]*x[i3,j2]*LL[i1,i2,i3]) = RHS[i1];