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: Mon, 20 Sep 2010 12:10:04 +0000 (UTC)
Organization: John D'Errico (1-3LEW5R)
Lines: 76
Message-ID: <i77ius$era$1@fred.mathworks.com>
References: <i735gj$782$1@fred.mathworks.com> <i73k78$k6d$1@fred.mathworks.com> <i75888$ais$1@fred.mathworks.com> <i75ch3$4in$1@fred.mathworks.com> <i775qs$rbo$1@fred.mathworks.com> <i77cb0$i6t$1@fred.mathworks.com> <i77fg8$6qg$1@fred.mathworks.com>
Reply-To: "John D'Errico" <woodchips@rochester.rr.com>
NNTP-Posting-Host: webapp-03-blr.mathworks.com
Content-Type: text/plain; charset=UTF-8; format=flowed
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1284984604 15210 172.30.248.38 (20 Sep 2010 12:10:04 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Mon, 20 Sep 2010 12:10:04 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 869215
Xref: news.mathworks.com comp.soft-sys.matlab:671497

"Alexander Pohl" <alexander.pohl@stfc.ac.uk> wrote in message <i77fg8$6qg$1@fred.mathworks.com>...
> "John D'Errico" <woodchips@rochester.rr.com> wrote in message <i77cb0$i6t$1@fred.mathworks.com>...
> > "Alexander Pohl" <alexander.pohl@stfc.ac.uk> wrote in message <i775qs$rbo$1@fred.mathworks.com>...
> > 
> > > Thank you very much for your very detailed instructions. As you already suspected I don't have the curve fitting toolbox with the fnint() routine. So would it be possible to send me your substitute?
> > > 
> > 
> > No. Actually, I specifically asked if you have the
> > SPLINES toolbox, which is where that function
> > would be found.
> > 
> > function ppint = myfnint(pp)
> > % integral of the piecewise cubic function pp, ppval(ppint,0) == 0
> > ppint = pp;
> > order = pp.order;
> > ppint.order = order + 1;
> > ppint.coefs = zeros(ppint.pieces, order + 1);
> > c = 0;
> > k = 1./(order:-1:1);
> > for i = 1:pp.pieces
> >   coefs = pp.coefs(i,:).*k;
> >   ppint.coefs(i,:) = [coefs,c];
> >   c = c + polyval([coefs,0],pp.breaks(i + 1) - pp.breaks(i));
> > end
> > 
> > This will work.
> > 
> > John
> 
> Dear John,
> 
> Thank you so much for your help. That code works beautifully now. One little note, the spline toolbox (with the fnint() function) is part of the curve fitting toolbox 3.0 which you can buy. I could not find the spline toolbox separately under the listed products.
> 

Ah yes. With the most recent release, they have moved it into
the curve fitting toolbox. I do need to download that release.
One more thing to do asap. Ugh.


> function [ P ] = integral( t ) 
> for i = 1:length(t)
>   P(i) = quadgk(@integrand, 0, t(i));
> end
>   function [ res ] = integrand( s )
>   res = exp(- s.^2 ) .* sin( s );
>   end
> end
> 
> function [ P ] = integral3( t ) 
> f = @(t) exp(-t.^2).*sin(t);
> fp = @(t) exp(-t.^2).*cos(t) - 2*t.*exp(-t.^2).*sin(t);
> slm = hermite2slm([t,f(t),fp(t)]);
> pp = slm2pp(slm);
> ppint = myfnint(pp);
> P = ppval(ppint,t);
> end
> 
> Benchmark:
> >> tbr = linspace(0,5,1000)';
> >> tic;P = integral(tbr)';toc
> Elapsed time is 1.142923 seconds.
> >> tic;P2 = integral3(tbr);toc
> Elapsed time is 0.061169 seconds.

In order to make this efficient, you don't want to 
recompute the splines every time you call the code.

For example, in the Fresnel sine/cosine integral
submission I just posted, I use a similar scheme
to generate a 7th degree Hermite interpolant. This
yields about 14 significant digits in the result. But
there I compute the splines ONCE. Then i use
persistent variables to keep the curve in memory
for when I will use it again.

John