Thread Subject:

How to optimise numerical integration?

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

"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

  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

"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

"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

"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

"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

"Alexander Pohl" <alexander.pohl@stfc.ac.uk> wrote in message <i74jk9$l34$1@fred.mathworks.com>...
> 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

"John D'Errico" <woodchips@rochester.rr.com> wrote in message <i75ch3$4in$1@fred.mathworks.com>...
> "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

Dear John,

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?

Best wishes,
Alexander

"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

"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.

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

"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

"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

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 );