Path: news.mathworks.com!not-for-mail
From: "John D'Errico" <woodchips@rochester.rr.com>
Newsgroups: comp.soft-sys.matlab
Subject: Re: How to optimise numerical integration?
Date: Sun, 19 Sep 2010 16:08:03 +0000 (UTC)
Organization: John D'Errico (1-3LEW5R)
Lines: 121
Message-ID: <i75ch3$4in$1@fred.mathworks.com>
References: <i735gj$782$1@fred.mathworks.com> <i73k78$k6d$1@fred.mathworks.com> <i75888$ais$1@fred.mathworks.com>
Reply-To: "John D'Errico" <woodchips@rochester.rr.com>
NNTP-Posting-Host: webapp-05-blr.mathworks.com
Content-Type: text/plain; charset=UTF-8; format=flowed
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1284912483 4695 172.30.248.35 (19 Sep 2010 16:08:03 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Sun, 19 Sep 2010 16:08:03 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 869215
Xref: news.mathworks.com comp.soft-sys.matlab:671266

"Alexander Pohl" <alexander.pohl@stfc.ac.uk> wrote in message <i75888$ais$1@fred.mathworks.com>...
> "John D'Errico" <woodchips@rochester.rr.com> wrote in message <i73k78$k6d$1@fred.mathworks.com>...
> > "Alexander Pohl" <alexander.pohl@stfc.ac.uk> wrote in message <i735gj$782$1@fred.mathworks.com>...
> > > Hi,
> > > 
> > > I would like to calculate the integral of exp(-t^2)sin(t) from 0 to t and make the following calculation as efficient as possible. Currently each point of the t vector has to be evaluated individually inside the 'for loop' which takes a long time.
> > > 
> > > Does anybody know how to vectorise the 'for loop' or any other method which would beat the current code in calculation time? I've got the optimisation toolbox and the NAG toolbox available if that's of any use.
> > >
> > 
> > They won't help much.
> > 
> > Personally, I'd be tempted to integrate a piecewise Hermite
> > interpolant. You know the derivatives of this function. They
> > are trivial to evaluate. So build a 5th or 7th order Hermite 
> > interpolant, then integrate that interpolant exactly.
> > 
> > John
> 
> Dear John,
> 
> If I approximate my function exp(-t^2)sin(t) with a Hermite polynomial or a spline function before integrating how do I integrate from 0 to t without the for loop? Performance wise this is currently much worse than my initial attempt.
> 
> function [ P ] = integral3( t ) 
> y = exp(- t.^2 ) .* sin( t );
> pp = spline(t, y);
> for i = 1:length(t)
>   P(i) = quad(@(t)ppval(pp,t), 0, t(i));
> end
> end
> 
> Could you give me a code example of how you would do it?
> 
> Many thanks for your help,
> Alexander

Wrong! Don't use a tool like quad to integrate a spline
function!

A cubic spline is just composed of cubic polynomial
segments. Cubic polynomials are rather trivial to
integrate.

By the way, spline does not produce a Hermite interpolant,
nor does pchip produce a useful hermite interpolant, at
least not one that will be accurate for you.

If I wanted to find an efficient scheme for a high accuracy
integration of a simple function like this, it would not use
quad at all. For example, try this...

The maximum value that we will integrate to...

>> tmax = 5;
>> tbr = linspace(0,tmax,100)';

I'll use a 5th order Hermite interpolant for this example.
This requires that I supply the function values at the
breaks, as well as the first and second derivatives. It
looks like I correctly got those derivatives off the back
of an envelope...

>> f = @(t) exp(-t.^2).*sin(t);
>> fp = @(t) exp(-t.^2).*cos(t) - 2*t.*exp(-t.^2).*sin(t);
>> fpp = @(t) -2*t.*exp(-t.^2).*cos(t) - exp(-t.^2).*sin(t) + ...
  -2*t.*fp(t) - 2*f(t);

Now, use a couple of tools from my SLM tool set,
the first one generates an interpolant from those
derivatives.

>> slm = hermite2slm([tbr,f(tbr),fp(tbr),fpp(tbr)]);

Next, convert it to a pp form.

>> pp = slm2pp(slm);

Integrate the interpolant, yielding a new function that
is quite efficient to evaluate at any upper limit. This
integration is easily done with fnint from the splines
toolbox, but the integration is trivial even without
that TB. (Ask if you don't have fnint, I can give you
a substitute.)

>> ppint = fnint(pp);
>> tic,int_pp = ppval(ppint,5);toc
Elapsed time is 0.000132 seconds.

>> tic,int_q = quadgk(f,0,5);toc
Elapsed time is 0.000811 seconds.

>> int_pp
int_pp =
          0.424436383489929

>> int_q
int_q =
          0.42443638350328

Note that the two agree to about 10 digits. Had I wanted
a more accurate approximation, I might have used a finer
set of knots, or a higher order interpolant. (I'm feeling too
lazy now to compute the third derivative.) With a few more
breaks in that spline, I get full double precision accuracy.

>> int_pp - int_q
ans =
     -1.33511535160835e-11

The real gain comes if I want to evaluate the integral at
many points.

>> tic,int_pp = ppval(ppint,0:.000001:5);toc
Elapsed time is 1.152340 seconds.

I updated hermite2slm and slm2pp recently to handle
higher order interpolants. I'll post the new version of
slm with those updates, so it will appear on Monday,
9/20/2010.

John