Newton Cotes: 1D

The idea is to make an integration exact on polynomials up to some order.  First recreating closed Newton-Cotes formulae

Evaluation of a polynomial at a bunch of points is multiplication of the coefficient list by the vandermonde matrix.

Evaluation of a polynomial at a a bunch of points is multiplication of the differentiated coeff list by the vm matrix. Easy to code.

The integration of a polynomial on the interval defined by xList (assuming it is ordered) is easy.  It is another matrix and then evaluation

Clear[vandermonde,xList,m]
vandermonde[xList_,m_]:=Module[
{i},
 Transpose[
Join[
{0*xList+1},
Table[(xList^i),{i,1,m-1}]
]
]
]
DMat[m_]:= SparseArray[
Table[{i,i+1}->i,{i,1,m-1}],
{m,m}]
IntMat[m_]:= SparseArray[
Table[{i+1,i}->(1/i),{i,1,m}],
{m+1,m}]

Clear[n,x,cList,VMat]
n=8;
xList=Range[0,1,1/2];
cList=RandomReal[{-1,1},n];
p[x_]=Sum[cList[[i]] (x^(i-1)),{i,1,n}];
P[x_]=Integrate[p[x],x];
VMat[n_]:=VMat[n]=vandermonde[xList,n];
{Norm[p[xList]-VMat[n].cList],
Norm[p'[xList]-VMat[n].DMat[n].cList],
Norm[p''[xList]-VMat[n].DMat[n].DMat[n].cList ],
Norm[p'''[xList]-VMat[n].MatrixPower[DMat[n],3].cList]
}
Norm[P[xList]-VMat[n+1].IntMat[n].cList]
Integrate[p[x],{x,xList[[1]],xList[[-1]]}]
(VMat[n+1][[-1]]-VMat[n+1][[1]]).IntMat[n].cList

So I now know how to do all the damm stuff analytically.  I can make sure I get the std formulae.

I need things which are exact on polynomials so for Newton Cotes All I do is make sure that some vector of weights dotted with the polynomial evaluated at the specifed points is the same as the analytic integral

First here are the Newton Cotes formulae.

Clear[w1,w2,w3]
xList=Range[0,1,1/2];
cList={c0,c1,c2,c3};
n=Length[cList];
SolveAlways[
Simplify[VMat[n].cList.{w1,w2,w3} == (VMat[n+1][[-1]]-VMat[n+1][[1]]).IntMat[n].cList],
{c0,c1,c2}]
SolveAlways[
Simplify[Transpose[VMat[n]].{w1,w2,w3}.cList == (VMat[n+1][[-1]]-VMat[n+1][[1]]).IntMat[n].cList],
{c0,c1,c2}]
Simplify[Transpose[VMat[n]].{w1,w2,w3}== (VMat[n+1][[-1]]-VMat[n+1][[1]]).IntMat[n]]
LinearSolve[Transpose[ VMat[n]], (VMat[n+1][[-1]]-VMat[n+1][[1]]).IntMat[n]]

Here is a table of NC formulae to show that it works.

Clear[VMat,n] 
Table[
Clear[VMat];
xList=Range[0,1,1/(n-1)];
VMat[n_]:=VMat[n]=vandermonde[xList,n];
NCFormulae=LinearSolve[Transpose[ VMat[n]], (VMat[n+1][[-1]]-VMat[n+1][[1]]).IntMat[n]],
{n,2,6}]

What would be interesting is a table of Newton-Cotes-Hermite formulae.

The first thing to do is make sure I have my head on straight.

Clear[VMat,n];
xList=Range[0,1,1/2];
cList={c0,c1,c2,c3};
n=Length[cList];
VMat[n_]:=VMat[n]=vandermonde[xList,n];
{w1,w2,w3,v1,v2,v3}/.SolveAlways[
Simplify[VMat[n].cList.{w1,w2,w3} + VMat[n].DMat[n].cList.{v1,v2,v3}== (VMat[n+1][[-1]]-VMat[n+1][[1]]).IntMat[n].cList],
cList][[1]]

Clear[VMat,n]
cList={c0,c1,c2,c3,c4};
n=Length[cList];
VMat[n_]:=VMat[n]=vandermonde[xList,n];
{w1,w2,w3,v1,v2,v3}/.SolveAlways[
Simplify[VMat[n].cList.{w1,w2,w3} + VMat[n].DMat[n].cList.{v1,v2,v3}== (VMat[n+1][[-1]]-VMat[n+1][[1]]).IntMat[n].cList],
cList][[1]]

Clear[VMat,n]
cList={c0,c1,c2,c3,c4,c5};
n=Length[cList];
VMat[n_]:=VMat[n]=vandermonde[xList,n];
{w1,w2,w3,v1,v2,v3}/.SolveAlways[
Simplify[VMat[n].cList.{w1,w2,w3} + VMat[n].DMat[n].cList.{v1,v2,v3}== (VMat[n+1][[-1]]-VMat[n+1][[1]]).IntMat[n].cList],
cList][[1]]

Now I need to build this as a matrix solve.

Table[
Clear[VMat];
xList=Range[0,1,1/2];
VMat[n_]:=VMat[n]=vandermonde[xList,n];
A=ArrayFlatten[{{Transpose[VMat[n]], Transpose[VMat[n].DMat[n]]}}];
b= (VMat[n+1][[-1]]-VMat[n+1][[1]]).IntMat[n];
{LinearSolve[A,b],NullSpace[A]},
{n,3,6}]

Table[
Clear[VMat];
xList=Range[0,1,1/3];
VMat[n_]:=VMat[n]=vandermonde[xList,n];
A=ArrayFlatten[{{Transpose[VMat[n]], Transpose[VMat[n].DMat[n]]}}];
b= (VMat[n+1][[-1]]-VMat[n+1][[1]]).IntMat[n];
{LinearSolve[A,b],NullSpace[A]},
{n,3,8}]

Including second derivatives.

Table[
Clear[VMat];
xList=Range[0,1,1/2];
VMat[n_]:=VMat[n]=vandermonde[xList,n];
A=ArrayFlatten[{{Transpose[VMat[n]], Transpose[VMat[n].DMat[n]],Transpose[VMat[n].DMat[n].DMat[n]]}}];
b= (VMat[n+1][[-1]]-VMat[n+1][[1]]).IntMat[n];
{LinearSolve[A,b],NullSpace[A]},
{n,3,9}]

Table[
Clear[VMat];
xList=Range[0,1,1/3];
VMat[n_]:=VMat[n]=vandermonde[xList,n];
A=ArrayFlatten[{{Transpose[VMat[n]], Transpose[VMat[n].DMat[n]],Transpose[VMat[n].DMat[n].DMat[n]]}}];
b= (VMat[n+1][[-1]]-VMat[n+1][[1]]).IntMat[n];
{LinearSolve[A,b],NullSpace[A]},
{n,3,12}]

Optimal Elements in Famillies

I have these families of integrators in 1D. There are several questions.  
1) Are they any good in comparison to the Newton-Cotes formulae?
2) How do they generalize to more than one dimension?

2D

I want to organize a polynomial basis in 2D.  I have two obvious choices: 
	A) bilinear, biquadratic, bicubic, biquartic etc.
	B) Complete Linear, Complete Quadratic, Complete Cubic etc.
I do not think the first one makes sense for triangles since I do not believe it is invariant under linear transformations.

Here is a 1 Dimensional Basis using symbol "x" and the easiest indexing version in Mathematica,

Clear[x,p,i]
Basis[x_,1,p_]:=Table[((x[1])^i),{i,0,p}]
p=3;
Temp1D=Basis[x,1,p]

One way to build a basis for the complete 2D quadratics is inductively.

Temp2D=Table[
Basis[x,1,p-i]((x[2])^i),
{i,0,p}]

Now it is easy to flatten the structure of this set but the structure might be usefull for building matrices indcutively.

Flatten[Temp2D]

When it is time to do it I can build 3D from 2D.

Temp1D.{c0,c1,c2,c3}
Temp2D[[1]].{c00,c01,c02,c03}+Temp2D[[2]].{c10,c11,c12}+Temp2D[[3]].{c20,c21}+Temp2D[[4]].{c30}
Flatten[Temp2D].Flatten[{{c00,c01,c02,c03},{c10,c11,c12},{c20,c21},{c30}}]

N Dimensions

I can easily write down a basis for polynomials in "n" dimensions.  I can even use restricted polynomial space appropriate for certain FEA analysis.

Basis[x_,1,p_]:=Table[((x[1])^i),{i,0,p}] 
Basis[x_,n_,p_]:= Sort[
Flatten[
Table[
Basis[x,n-1,i]((x[n])^(p-i)),
{i,0,p}]
],
Exponent[#1/.x[i_]->xx,xx]>Exponent[#2/.x[i_]->xx,xx]&
]
Basis[x,1,2].Array[c,3]

Basis[x_,1,p_]:=Table[((x[1])^i),{i,0,p}] 
Basis[x_,n_,p_]:= Table[
Basis[x,n-1,i]((x[n])^(p-i)),
{i,p,0,-1}]
Basis[x,1,2].Array[c,3]
D[Basis[x,1,2],x[1]].Array[c,3]
D[Basis[x,2,2],x[2]]

D[Basis[x,2,5],{Array[x,2]}]

Clear[p]
Product[ ((x[k])^(p[k])),{k,1,3}]

Basis[x]

Array[y,3]

Table[
IntegerDigits[i,3,3],
{i,0,12}]

IntegerPartitions[5,{3}]

?IntegerPartitions