% matlab script to derive the 2-point Gauss-Laguerre quadrature rule % and use it on an example % inelegantly set up equations from method of undetermined coefficients % and solve clear all close all format long eq1 = 'w1*1 + w2*1 = 1'; eq2 = 'w1*x1 + w2*x2 = 1'; eq3 = 'w1*x1^2 + w2*x2^2 = 2'; eq4 = 'w1*x1^3 + w2*x2^3 = 6'; [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 = @(x) exp(x).*log(1+exp(-x)); % use Gauss-Laguerre quadrature to approximate it glq = w1*feval(f,x1) + w2*feval(f,x2) % now let's find the exact answer and see what the error is ... syms x I = double(int(log(1+exp(-x)),0,Inf)) error = I - glq pause % or we can estimate the neglected tail % trying integral(f,0,Inf) or integral(f,0,realmax) % won't work! f = @(x) log(1+exp(-x)); Q10 = integral(f,0,10) Q20 = integral(f,0,20) Q30 = integral(f,0,30) % from here we see that we have convergence to 8 decimal places % we can also perform a change of variable to make the interval finite F = @(t) log(1+t)./t; integral(F,0,1) % check for problems at t=0 ezplot(F,0,1) integral(F,realmin,1)