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