Re: definite numerical integration

[David Bateman]

roberto wrote:

> On 9/21/05, David Bateman <address@hidden> wrote:
>  
>
> > roberto wrote:
> >    
> >
> > > 3. doing this, evaluating the integral by varying r* in [0,R] (e.g. R = 110)
> > > i'll have a function F(r*,a) where a should be obtained by imposing
> > > some boundary condition;
> > >
> > > i looked for some help into Octave like
> > > http://octave.sourceforge.net/index/Q.html
> > >
> > > but in libraries like:
> > > [v, ier, nfun, err] = quad (f, a, b, tol, sing)
> > >
> > > the function to be integrated are to be inserted by their explicit formulation
> > > which i don't have now, as already stated above;
> > >
> > >
> > >
> > >      
> > I don't quite see the problem as in 1) you saw you have a function to
> > numerically calculate the values in a range [0,r_max] and then in 3) you
> > state you can use a explicit formulation
> >
> >    
> sorry, but actually i do not have an explicit formulation of f(r) like 
> function y = fun (x)
> # Nasty integrable singularity to be unkind to quad and show octave's
>   better here :-)
>     y  = 1/x;
> endfunction
>
> but instead i have only values of f(r) in a finite set of points r,
> maybe this is useful to solve the problem
>
>
>
>  
Then all you can do is something like

x0 = ??;
y0 =  f(x0);
a = ??;
b = ??;

# x0 must be montonically increasing
[x0, idx] = sort(x0);
y(idx) = y0;

# Treat core of the integral with simpson's rule
Istart = find (x0 >= a)(1);
Istop = find (x0 <= b)(end);
Int = sum(0.5 * (y0((Istart+1):Istop) + y0(Istart:(Istop-1))) .* 
(x0((Istart+1):Istop) - x0(Istart:(Istop-1))));

# Treat the tail
if (x0(Istart) != a)
 Int += 0.5 * (x0(Istart) - a) * (y0(Istart) + y0(Istart-1));
endif
if (x0(Istop) != b)
 Int += 0.5 * (b - x0(Istop)) * (y0(Istop+1) + y0(Istop));
endif

Though you can't get the estimate of the error in the integral in this 
manner, but you know its pretty much the best you can do, unless the 
values of x0 might be chosen in such a way as to fall on the absicca of 
a Guassian quadrature, in which case a guassian rule might be used to 
replace simpson's rule.

Testing the above with

x0=0:0.0333:1;
y0=x0;
a=0
b=0.5;

Should give 1/8... It gives "Int =  1.250082e-01", which is pretty good.

D.

--
David Bateman                                address@hidden
Motorola Labs - Paris                        +33 1 69 35 48 04 (Ph) 
Parc Les Algorithmes, Commune de St Aubin    +33 1 69 35 77 01 (Fax) 
91193 Gif-Sur-Yvette FRANCE

The information contained in this communication has been classified as: 

[x] General Business Information 
[ ] Motorola Internal Use Only 
[ ] Motorola Confidential Proprietary



-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------