On 2 Jun 2009, at 22:15, Carlo de Falco wrote: > > On 2 Jun 2009, at 21:12, dmelliott wrote: > >> >> Dear Søren, >> >> More thoughts: >> >> >> The idea of this routine is that it integrates existing data in >> its most >> ususal form. That is, it does not need a function to call, and, >> since, for >> better or for worse, most data comes on a domain of equally sized >> steps, >> this is what it uses. This is an attempt to make it compatable >> whith most >> experimantal data, and the outputs of other routines, e.g. sales by >> week, >> Fourie analysed amplitudes by frequency, actuarial data by age, >> etc. The >> patrial interval integrals, including cumulative, can be achieved >> by looping >> or otherwise manipulating the endpoint values. >> >> The best way to see its worth is to compare it to the available >> routine >> for existing data, by editing the "integrator_verifier" to use the >> "trapz" >> routine. The latter is only of first order, and even for as many >> as 500 >> points on [0,1] yields a 2.01e-6 relative error for the next higher >> order, >> the second. For higher orders higher than second it is even more >> useless. >> There probably should a caution appended to the instructions about >> this. > > Douglas, > > It seems to me there is one detail you are missing here: > experimental data are not always polynomials, actually most often > you don't even have a way to tell wether they are smooth, if all you > have are equispaced samples. > > Using a high order quadrature rule only makes sense if you know a- > priori your data can be > accurately approximated by a polynomial, e.g. by Taylor expansion. > If you have no a-priori knowledge about the smoothness of your data, > using a low order quadrature rule is a safer bet. For example > consider this well known example: > > >> x = linspace (-1, 1, 10); > >> f = @(x) 1./(1+25.*x.^2); > >> F = f(x); > >> integrator (F, -1, 1) > ans = 0.47972 > >> sum ((F(2:end)+F(1:end-1)).*diff(x)/2) > ans = 0.54437 > > while the exact value of the integral is > > >> .4 * atan (5) > ans = 0.54936 > >> quad(f, -1, 1) > ans = 0.54936 > > so in this case even the simple trapezoidal rule beats your 10th > order method because it's based on approximating f by means of a > highly-order and thus highly oscillating polynomial. > > So I believe your function would actually be more useful if the > order of the quadrature rule to use was to be given as an input > option rather than be set automatically depending on the size of the > sample. > > HTH, > c.
Another, even simpler example : >> x = linspace (-1, 1, 15); >> f = @(x) abs(x); >> F = f(x); >> integrator (F, -1, 1) ans = 9.92164981733434e-01 >> sum ((F(2:end)+F(1:end-1)).*diff(x)/2) ans = 1.00000000000000e+00 >> quad(@(x) abs(x), -1, 1) ans = 1.00000000000000e+00 >> In this case, quadrature by the trapezoidal rule is exact as the integrand function is piece-wise linear while a higher order rule suffers because the first derivative has a jump so the 2nd derivative is unbounded. c. ------------------------------------------------------------------------------ OpenSolaris 2009.06 is a cutting edge operating system for enterprises looking to deploy the next generation of Solaris that includes the latest innovations from Sun and the OpenSource community. Download a copy and enjoy capabilities such as Networking, Storage and Virtualization. Go to: http://p.sf.net/sfu/opensolaris-get _______________________________________________ Octave-dev mailing list [email protected] https://lists.sourceforge.net/lists/listinfo/octave-dev
