Re: Better quadrature routine in octave

[Pedro Gonnet]

Hello again,

Ok, thanks for the tip regarding "x". I've got everything working fine
now except for the treatment of non-finite intervals, since this
requires me to wrap the integrand in another function.

In Octave, the code that did this was pretty straight-forward:

        f = @(x) f ( tan ( pi / 2 * x ) ) .* ( 1 + tan ( pi * x / 2 ).^2 ) * pi 
/ 2;

In the Oct-File, the integrand "f" is a octave_function pointer which I
extracted from the arguments as per section A.1.8 of the manual. I tried
creating a wrapper following the example at the end of that section and
couldn't get it to work since the example is buggy (std::octave?) and I
couldn't find any documentation for "extract_function". 

I don't even know if this is the best way to implement the above wrapper
in an Oct-File. Any ideas how this should be done? 

This is the only thing missing for the Oct-File version of cquad --
everything else works and is ready for prime-time.

Cheers and thanks,
Pedro



On Sat, 2010-04-24 at 17:14 -0400, John W. Eaton wrote:
> On 24-Apr-2010, Pedro Gonnet wrote:
> 
> | I'm currently re-writing the C-language version as an Oct-file. The
> | documentation is a bit shaky, but I can get by on the examples except
> | for one thing: When I call the integrand (I'm using an octave_function
> | *fcn), I need to generate a vector of values "x" to compute "y = f(x)".
> | I have to do this often and for vectors "x" of length 5, 17, 33, 3, 7,
> | 15 and 31.
> | 
> | I'm guessing the most efficient way would be do declare one Matrix
> | object of length 33 and then, whenever I call the integrand, fill it
> | with the values of "x" I need and somehow tell the Matrix to truncate
> | its own length.
> 
> If you resize a matrix, you may end up generating a copy.  So I would
> just create a matrix of the size you need each time you need it.
> 
> Also, how are you accessing the pointer to the octave_function object?
> It will probaby be more robust to just call feval.
> 
> If you want help, maybe it would be best to post what you have so far
> so someone could comment on it.  The help list would be a better place
> for that kind of question.
> 
> jwe