Re: Numerical integration: Simpson and trapezoidal method

[c.]

On 4 Oct 2010, at 14:27, David Lucas wrote:

> Hello,
>
> i've been looking for how to perform a numerical integration with  
> the Simpson's Rule and i found it:
>
> e.g.:
>         function y = f(x)
>                 y = sqrt(sin(x))
>         endfunction
>         
>         quadv("f", 1, pi/2)
>
> It gives me the result of the integration of sqrt(sin(x)) from 1 to  
> PI/2.
>
> The problem is that i also want to set the number of sub-intervals.

quadv adaptively computes the number of sub-intervals automatically in  
order to reduce the integration,
so you cannot set the number of subdivisions a-priori...

I don't know if there is a function to approximate an integral with  
the composite Simpson rule
for a given number of subdivisions but you can easily do it yourself  
as follows:

f = @(x) sqrt (sin (x)); a = 1; b = pi/2; N = 4;
x = linspace (a, b, N+1); x1 = x(1:end-1); x3 = x(2:end); x2 = (x1+x3)/ 
2;
q = sum ((b - a) / ( N * 6) * (f(x1) + 4 * f(x2) + f(x3)))

> In Maple i would do it like this(4 sub-intervals/partitions):
>
> 	ApproximateInt(sqrt(sin(x)), x = 1 .. 3.141592654*(1/2), method =  
> simpson, partition = 4)
>
> What about trapezoidal method? Does trapz have the same syntax as  
> quadv?

no, for trapz you pass directly the list of points where your function  
is to be evaluated and the corresponding values of f(x)

f = @(x) sqrt (sin (x)); a = 1; b = pi/2; N = 4;
x = linspace (a, b, N+1);
trapz (x, f(x))

> Also, is there any way to plot something like Maple does with the  
> following sentence?:
>
> ApproximateInt(sqrt(sin(x)), x = 1 .. 3.141592654*(1/2), method =  
> simpson, partition = 4, output = plot)
>
> Result: http://img148.imageshack.us/img148/2031/dibujoqw.jpg


plot (x, f(x))
hold on
stem (x, f(x), 'r')

> Thank you.

HTH,
c.