# Thread Subject: How to optimise numerical integration?

 Subject: How to optimise numerical integration? From: Alexander Pohl Date: 18 Sep, 2010 19:56:04 Message: 1 of 13 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. Here is the code: t = 0:0.001:5 tic; P = integral(t); toc; Elapsed time is 3.766351 seconds. (2 GHz intel) integral.m file: function [ P ] = integral( t ) for i = 1:length(t) %<---- problem%   P(i) = quad(@integrand, 0, t(i)); end   function [ res ] = integrand( s )   res = exp(- s.^2 ) .* sin( s );   end end Thanks for any help, Alexander
 Subject: How to optimise numerical integration? From: Roger Stafford Date: 18 Sep, 2010 21:31:03 Message: 2 of 13 "Alexander Pohl" wrote in message ... > 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. > > Here is the code: > t = 0:0.001:5 > tic; P = integral(t); toc; > Elapsed time is 3.766351 seconds. (2 GHz intel) > > integral.m file: > function [ P ] = integral( t ) > for i = 1:length(t) %<---- problem% > P(i) = quad(@integrand, 0, t(i)); > end > function [ res ] = integrand( s ) > res = exp(- s.^2 ) .* sin( s ); > end > end > > Thanks for any help, > Alexander - - - - - - - - -   What you have here is an inherently inefficient procedure. Five thousand and one times you are starting the integration back at 0 and up to an increasing value of t. If each time you integrated only from t(i) to t(i+1), the elapsed time should be far less. Then to get P you could do a simple cumsum on the individual integral values you had found. Roger Stafford
 Subject: How to optimise numerical integration? From: Roger Stafford Date: 18 Sep, 2010 22:06:03 Message: 3 of 13 I should point out that with points as closely-spaced as you have here (0:0.001:5) you ought to obtain a very decent accuracy using cumtrapz rather than quad, and that should make it even faster. Roger Stafford
 Subject: How to optimise numerical integration? From: John D'Errico Date: 19 Sep, 2010 00:07:04 Message: 4 of 13 "Alexander Pohl" wrote in message ... > 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
 Subject: How to optimise numerical integration? From: Alexander Pohl Date: 19 Sep, 2010 09:03:05 Message: 5 of 13 "Roger Stafford" wrote in message ... > "Alexander Pohl" wrote in message ... > > 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. > > > > Here is the code: > > t = 0:0.001:5 > > tic; P = integral(t); toc; > > Elapsed time is 3.766351 seconds. (2 GHz intel) > > > > integral.m file: > > function [ P ] = integral( t ) > > for i = 1:length(t) %<---- problem% > > P(i) = quad(@integrand, 0, t(i)); > > end > > function [ res ] = integrand( s ) > > res = exp(- s.^2 ) .* sin( s ); > > end > > end > > > > Thanks for any help, > > Alexander > - - - - - - - - - > What you have here is an inherently inefficient procedure. Five thousand and one times you are starting the integration back at 0 and up to an increasing value of t. If each time you integrated only from t(i) to t(i+1), the elapsed time should be far less. Then to get P you could do a simple cumsum on the individual integral values you had found. > > Roger Stafford Dear Roger, I've tried to implement your suggestion: function [ Y ] = integral2( t ) for i = 1:length(t)   if (i == length(t))       P(i) = 0;   else       P(i) = quad(@integrand, t(i), t(i+1));   end end Y = cumsum(P);   function [ res ] = integrand( s )   res = exp(- s.^2 ) .* sin( s );   end end I've used the actual t vector of my data to test: t = 0.104:0.016:31.720; P1 = integral(t) P2 = integral2(t) The two vectors P1 and P2 are not identical. Could you give me a code example how to implement the piecewise integration? Alexander
 Subject: How to optimise numerical integration? From: Alexander Pohl Date: 19 Sep, 2010 14:55:04 Message: 6 of 13 "John D'Errico" wrote in message ... > "Alexander Pohl" wrote in message ... > > 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
 Subject: How to optimise numerical integration? From: John D'Errico Date: 19 Sep, 2010 16:08:03 Message: 7 of 13 "Alexander Pohl" wrote in message ... > "John D'Errico" wrote in message ... > > "Alexander Pohl" wrote in message ... > > > 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
 Subject: How to optimise numerical integration? From: Roger Stafford Date: 19 Sep, 2010 18:51:03 Message: 8 of 13 "Alexander Pohl" wrote in message ... > Dear Roger, > > I've tried to implement your suggestion: > > function [ Y ] = integral2( t ) > for i = 1:length(t) > if (i == length(t)) > P(i) = 0; > else > P(i) = quad(@integrand, t(i), t(i+1)); > end > end > Y = cumsum(P); > function [ res ] = integrand( s ) > res = exp(- s.^2 ) .* sin( s ); > end > end > > I've used the actual t vector of my data to test: > t = 0.104:0.016:31.720; > P1 = integral(t) > P2 = integral2(t) > > The two vectors P1 and P2 are not identical. Could you give me a code example how to implement the piecewise integration? > > Alexander - - - - - - - - -   If you have done "integral" as you did originally:  P(i) = quad(@integrand, 0, t(i)) integrating from 0 to t(i), then they couldn't be expected to match. The first P(1) is the integral from 0 to 0.104, whereas the first integral in "integral2" is the integral from 0.104 to 0.104+0.016 . It never computes the integral from 0 to 0.104 .   In any event it is unrealistic to expect identical results. Round-off errors will inevitably produce some differences and 'quad' makes errors of its own. What is important here is whether there was any appreciable improvement in execution time. As I said before integrating each time starting at zero up to gradually increasing values of t is inherently inefficient.   You should take John's advice seriously, Alexander. He has explained a method that allows much more accuracy in integrating between given successive steps in t than with trapezoidal integration, or it allows many fewer steps to be taken to achieve the same accuracy. The same is also true with higher order type integration schemes that approximate with cubic or fifth order polynomials, but utilizing the known derivatives of the function using Hermite interpolants makes things even more efficient. As they are at present, the matlab quadrature functions, quad, quadgk, etc. are incapable of utilizing the derivatives of the integrand functions they are given and must compensate by taking larger numbers of samples of those functions. Roger Stafford
 Subject: How to optimise numerical integration? From: John D'Errico Date: 20 Sep, 2010 10:17:04 Message: 10 of 13 "Alexander Pohl" wrote in message ... > 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
 Subject: How to optimise numerical integration? From: Alexander Pohl Date: 20 Sep, 2010 11:11:04 Message: 11 of 13 "John D'Errico" wrote in message ... > "Alexander Pohl" wrote in message ... > > > 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. 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. Best wishes, Alexander
 Subject: How to optimise numerical integration? From: John D'Errico Date: 20 Sep, 2010 12:10:04 Message: 12 of 13 "Alexander Pohl" wrote in message ... > "John D'Errico" wrote in message ... > > "Alexander Pohl" wrote in message ... > > > > > 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
 Subject: How to optimise numerical integration? From: Jonas Lundgren Date: 20 Sep, 2010 13:37:04 Message: 13 of 13 "Alexander Pohl" wrote in message ... > 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. > > Here is the code: > t = 0:0.001:5 > tic; P = integral(t); toc; > Elapsed time is 3.766351 seconds. (2 GHz intel) > > integral.m file: > function [ P ] = integral( t ) > for i = 1:length(t) %<---- problem% > P(i) = quad(@integrand, 0, t(i)); > end > function [ res ] = integrand( s ) > res = exp(- s.^2 ) .* sin( s ); > end > end > > Thanks for any help, > Alexander Hi Alexander, A fast and accurate alternative using cumsum (as suggested by Roger Stafford) with a Gauss-Legendre quadrature (see also submission #4540 on the File Exchange). /Jonas function [ P ] = integral( t ) % 3-point Gauss-Legendre quadrature xq = [-sqrt(3/5),0,sqrt(3/5)]; wq = [5/9;8/9;5/9]; % Quadrature points hh = (t(2)-t(1))/2; tmid = t(1:end-1) + hh; tquad = bsxfun(@plus,tmid(:),hh*xq); % Integral P = integrand(tquad)*(wq*hh); P = [0; cumsum(P)]; function [ res ] = integrand( s ) res = exp(- s.^2 ) .* sin( s );

### Everyone's Tags:

Separated by commas
Ex.: root locus, bode

### What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.