Path: news.mathworks.com!not-for-mail From: <HIDDEN> Newsgroups: comp.soft-sys.matlab Subject: Re: How to optimise numerical integration? Date: Sun, 19 Sep 2010 09:03:05 +0000 (UTC) Organization: Rutheford Appleton Lab Lines: 56 Message-ID: <i74jk9$l34$1@fred.mathworks.com> References: <i735gj$782$1@fred.mathworks.com> <i73b2n$k58$1@fred.mathworks.com> Reply-To: <HIDDEN> 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 1284886985 21604 172.30.248.38 (19 Sep 2010 09:03:05 GMT) X-Complaints-To: news@mathworks.com NNTP-Posting-Date: Sun, 19 Sep 2010 09:03:05 +0000 (UTC) X-Newsreader: MATLAB Central Newsreader 2508071 Xref: news.mathworks.com comp.soft-sys.matlab:671218 "Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <i73b2n$k58$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. > > > > 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