[Math] Numerical integration in Matlab (Gaussian 3 point quadrature)

integrationMATLABnumerical methods

Write a Matlab function that applies the Gauss three point rule to N sub-intervals of $[a, b].$ The input parameters should be the name of the function being integrated, $a, b,$ and $N$.

Attempt:

For $\int^b_a f(x) dx$ the Gaussian 3 point quadrature is given by:

$$\frac{b-a}{18} \left( 5f \left( \frac{a+b}{2}- \sqrt{\frac{3}{5}} \frac{b-a}{2} \right) +8f \left( \frac{a+b}{2} \right) + 5 f \left( \frac{a+b}{2}+ \sqrt{\frac{3}{5}} \frac{b-a}{2} \right) \right).$$

Here is my Matlab code that uses this equation to approximate $\int^b_a f(x) dx$:

function t = gauss3(func,a,b)
j=(((a+b)/2)-sqrt(3/5)*((b-a)/2));
k=((a+b)/2);
l=(((a+b)/2)+sqrt(3/5)*((b-a)/2));
t=((b-a)/18)*(5*feval(func,j) + 8*feval(func,k) + 5*feval(func,l))

function y=f(x)
y=exp(-x^2);
> t = gauss3(@f,0,1)

So how can I modify this code to take into account the sub-intervals $N$? The function should start as function t = gauss3(func,a,b,N).

P.S. In the code above I am trying to approximate the integral $\int^1_0 e^{-x^2} \ dx,$ using the Gaussian three point rule with error $<10^{-10}.$

Edit:

function y = gaussq(func,a,b,N);

h = (b-a)/N;
y = 0;
for i = 1:N,
    y = y+gauss3(func,a+(i-1)*h,a+i*h);
end
% y0 = 100;
% for N = logspace(0,6,7),
%     y1 = gaussq(@f,a,b,N);
%     err = abs((y1-y0)/y0);
%     if err < 1.0e-10,
%         break;
%     end
%     y0 = y1;
% end

Where:

function y=f(x)
y=exp(-x^2);

Best Answer

I wouldn't modify the code except to suppress output on the last line of gauss3.m. All you need is a function that chops up $[a,b]$ into $N$ subintervals and invokes gauss3.m $N$ times and adds up the results:

% gaussq.m

function y = gaussq(func,a,b,N);

h = (b-a)/N;
y = 0;
for i = 1:N,
    y = y+gauss3(func,a+(i-1)*h,a+i*h);
end

Then you want a script that invokes the above for various values of $N$ and compares results until it's satisfied that the error criterion is met:

% testg.m

a = 0;
b = 1;
y0 = 100;
for N = logspace(0,6,7),
    y1 = gaussq(@f,a,b,N);
    err = abs((y1-y0)/y0);
    if err < 1.0e-10,
        break;
    end
    y0 = y1;
end
fprintf('Integral = %17.15f, intervals = %d, error = %e\n',y1,N,err);

Hopefully after you analyze the above example to see how it operates you will be able to write your own, much better, code.

Related Question