kill(all);

integro in forma esatta un generico polinomio di grado 3, 
dopo averlo definito

poli(pluto) := a*pluto^3 + b*pluto^2 + c*pluto + d;

INTpoli : integrate(poli(xi),xi,-1,1);

ricavo l'integrale esatto

INTpoli : fullratsimp ( INTpoli );

verifico che i termini dispari diano contributo nullo, essendo l'intervallo di integrazione simmetrico.

pongo l'integrale esatto uguale alla stima per campionamento
a due punti; NON suppongo simmetrica la disposizione dei punti di campionamento

residuo : INTpoli - ( 

  w1*poli(xi1) 
+ w2*poli(xi2) 

);

imposto un sistema di equazioni *NON*lineari per ricavare 
le incognite w1,w2,xi1,xi2

eqs :
[
  diff(residuo, a)=0,
  diff(residuo, b)=0,
  diff(residuo, c)=0,
  diff(residuo, d)=0
];

risolvo il sistema di equazioni NONlineari

sols : solve( eqs , [w1,w2,xi1,xi2] ) ;

sols[1];

sols[2];

definisco le incognite sulla base di una delle due soluzioni

w1  : ev( w1,sols[1]);
w2  : ev( w2,sols[1]);
xi1 : ev(xi1,sols[1]);
xi2 : ev(xi2,sols[1]);

controllo che tali coefficienti, oltre a rendere il 
residuo costante in a,b,c,d, lo rendano anche nullo

ev( residuo );

fullratsimp(%);

provo ad estendere la formulazione ad intervalli generici, a <= x <= b

nel caso specifico provo ad integrare sin(x) con 0 <= x <= %pi

mappatura x(xi)

x(xi) := a*xi+b;

linsolve([ x(-1)=0 , x(1)=%pi ] , [a,b] ), globalsolve=true;

x(xi);

calcolo dx(xi)/dxi

define(
    dxdxi(xi),
    diff(x(xi),xi)
);

definisco l'integrando come funzione di xi

define( 
        integrando(xi) , 
        sin(x(xi)) * dxdxi(xi)
);

calcolo l'integrale secondo quadratura gaussiana e in forma esatta

quadgauss : w1 * integrando(xi1) + w2 * integrando(xi2);

esatto:     integrate(integrando(xi),xi,-1,1);

verifico l'errore

[
    esatto,
    quadgauss,
    (esatto-quadgauss)/esatto
]
,numer;