% matlab script to demonstrate use of Gauss quadrature clear all close all % first derive the 2-point Gauss quadrature rule eq1 = 'w1*1 + w2*1 = 2'; eq2 = 'x1*w1 + x2*w2 = 0'; eq3 = 'x1^2*w1 + x2^2*w2 = 2/3'; eq4 = 'x1^3*w1 + x2^3*w2 = 0'; [w1,w2,x1,x2] = solve(eq1,eq2,eq3,eq4) pause % note: there are two solutions % we pick the one where x1 < x2 [x1,index] = min(double(x1)) w1 = double(w1(index)) x2 = double(x2(index)) w2 = double(w2(index)) % define the integrand f = @(t) exp(-((t+1)./2).^2) % use Gauss-Legendre quadrature to approximate it gq = (w1*feval(f,x1) + w2*feval(f,x2)); % don't forget the scaling! gq = gq/2 % now let's find the exact answer and see what the error is ... syms x I = double(int(exp(-x^2),0,1)) errorGQ = I - gq % compare to Simpson's rule a = -1; b = 1; c = (a+b)/2; s = (b-a)/6*(feval(f,a)+4*feval(f,c)+feval(f,b)); % don't forget the scaling! s = s/2 errorS = I - s