MATLAB Central Newsreader - trapz, numerical integration

trapz, numerical integration

leo nidas — Tue, 02 Jun 2009 22:27:02 -0400

Hi there,

I have a code in a for loop. At the end of each iteration a (matematical) function, let
f(x), comes up which is always decreasing. x>0 and f(x)->0 as x->Inf, also f(x)>0. Actually it is going exponentially to zero.

My goal is to integrate numerically (using the trapz function) the function from 0 to Inf.
Obviously I cannot insert Inf into the trapz (at least I think so, because it will at least take forever :-) ).

The thing is that I do not know a priori the exact form of the function to be integrated so as to insert a logical "large" number as an upper bound in the trapz. Till now I visualized some repetitions of the function that came up and the value 60 seemed more that large. But I want to be 100% sure that none of my functions gets away without proper integration.

Is there any trick or any other function that "detects" the logical point up to a function should be integrated, if one wants to integrate it till Inf?

Should I create a condition such as if f(x(i))<0.000001 then x(i) is my number? Or is there any better way?

Thanx in advance!!

P.S. Anallytical integration is not possible.

Re: trapz, numerical integration

John D'Errico — Tue, 02 Jun 2009 23:15:04 -0400

"leo nidas" <bleonidas25@yahoo.gr> wrote in message <h048vl$2v$1@fred.mathworks.com>...
>
> Hi there,
>
> I have a code in a for loop. At the end of each iteration a (matematical) function, let
> f(x), comes up which is always decreasing. x>0 and f(x)->0 as x->Inf, also f(x)>0. Actually it is going exponentially to zero.
>
> My goal is to integrate numerically (using the trapz function) the function from 0 to Inf.
> Obviously I cannot insert Inf into the trapz (at least I think so, because it will at least take forever :-) ).
>
> The thing is that I do not know a priori the exact form of the function to be integrated so as to insert a logical "large" number as an upper bound in the trapz. Till now I visualized some repetitions of the function that came up and the value 60 seemed more that large. But I want to be 100% sure that none of my functions gets away without proper integration.
>
> Is there any trick or any other function that "detects" the logical point up to a function should be integrated, if one wants to integrate it till Inf?
>
> Should I create a condition such as if f(x(i))<0.000001 then x(i) is my number? Or is there any better way?
>
>
> Thanx in advance!!
>
> P.S. Anallytical integration is not possible.

Consider a Gauss-Laguerre or generalized Gauss-Laguerre
quadrature. You can surely find something on the FEX.

Or, transform your problem, using a mapping to a
finite interval.

Then use your favorite quadrature tool.

John

Re: trapz, numerical integration

leo nidas — Sun, 07 Jun 2009 08:40:17 -0400

Thanx John,

I downloaded quadgr that seemed better than others if one wanted to integrate till Inf.

But consider the case:

b0=-5;
b1=3;
m=3

pc= quadgr(@(t) exp(b0+b1.*t)./(1+exp(b0+b1.*t)).*1/m.*exp(-t./m),0,Inf)

%plot the function I want to integrate:
gr=0:0.01:30;
plot(gr,exp(b0+b1.*gr)./(1+exp(b0+b1.*gr)).*1/m.*exp(-gr./m))

Quadgr doesn't seem to work here. I know that if I use trapz I could get a proper approximation but I need the integration untill Inf because my code is in a for loop and a new function (of the above form) to be integrated comes up. I.e. tha b0, b1, m may be different in each iteration with m>0.

Thanx again for any help.

Re: trapz, numerical integration

John D'Errico — Sun, 07 Jun 2009 12:05:02 -0400

"leo nidas" <bleonidas25@yahoo.gr> wrote in message <h0fudh$pfo$1@fred.mathworks.com>...
>
> Thanx John,
>
> I downloaded quadgr that seemed better than others if one wanted to integrate till Inf.
>
> But consider the case:
>
> b0=-5;
> b1=3;
> m=3
>
> pc= quadgr(@(t) exp(b0+b1.*t)./(1+exp(b0+b1.*t)).*1/m.*exp(-t./m),0,Inf)
>
> %plot the function I want to integrate:
> gr=0:0.01:30;
> plot(gr,exp(b0+b1.*gr)./(1+exp(b0+b1.*gr)).*1/m.*exp(-gr./m))
>
>
>
> Quadgr doesn't seem to work here. I know that if I use trapz I could get a proper approximation but I need the integration untill Inf because my code is in a for loop and a new function (of the above form) to be integrated comes up. I.e. tha b0, b1, m may be different in each iteration with m>0.
>
> Thanx again for any help.

I had originally suggested the use of a Gauss-Laguerre
numerical integration. If your function is this:

fun = @(t) exp(b0+b1.*t)./(1+exp(b0+b1.*t)).*1/m.*exp(-t./m)

then it can be broken down into two pieces.

    exp(b0+b1.*t)./(1+exp(b0+b1.*t)).*1/m

and

    exp(-t./m)

I'll look at the part which is not exponential to see if a
Gauss-Laguerre will work well.

K =@(t) exp(b0+b1.*t)./(1+exp(b0+b1.*t))./m;
ezplot(K)

It is not really that polynomial, although it is well behaved,
with a sigmoidal shape.

I then tried a Gauss-Laguerre quadrature on this function. It
is not polynomial enough that a low order Gauss-Laguerre
works, and roots fails above 14 nodes to generate any higher
order Gauss-Laguerre rules than that. (I should have used a
better routine to solve for the Gauss rule. No matter, since
there are better choices to be found.) Look instead at a
numerical tool like quadgk.

quadgk(fun,0,30,'reltol',1e-12)
ans =
         0.584687911313577

quadgk(fun,0,100,'reltol',1e-12)
ans =
         0.584733311243336

quadgk(fun,0,200,'reltol',1e-12)
ans =
         0.584733311243339

Or, try quadl. Note that quadl used more than 1200
function evaluations to get this result, but I had set
the tolerance rather small.

[I,fev] = quadl(fun,0,200,1e-12)
I =
         0.584733311243339
fev =
        1248

The default tolerance for quadl gave me a reasonable
result at a cost of only 168 function evaluations. In fact
this is far better accuracy than I would get from trapz
anyway.

[I,fev] = quadl(fun,0,200)
I =
         0.584733310572222
fev =
   168

How well does trapz do? At a cost of 3000 function
evaluations, you only get about 4 digits of accuracy,
whereas the tools like quadl or quadgk will yield
twice as many correct digits for vastly less work.

t = 0:0.01:30;
trapz(t,fun(t))
ans =
         0.584687862068707

So I'm not at all sure why you have chosen to use
trapz here, perhaps only that you did not know these
other tools exist in MATLAB. Trapezoidal integration
is rarely a good choice.

Can we do better yet? Certainly so, if I knew that your
function would always have the same essential character
of a sigmoid multiplied by a negative exponential.

I'll stop here for now, unless you wanted me to explore
this last avenue more deeply. Ask again if so.

John

Re: trapz, numerical integration

leo nidas — Sun, 07 Jun 2009 16:53:01 -0400

> I had originally suggested the use of a Gauss-Laguerre
> numerical integration. If your function is this:
>
> fun = @(t) exp(b0+b1.*t)./(1+exp(b0+b1.*t)).*1/m.*exp(-t./m)
>
> then it can be broken down into two pieces.
>
> exp(b0+b1.*t)./(1+exp(b0+b1.*t)).*1/m
>
> and
>
> exp(-t./m)
>
> I'll look at the part which is not exponential to see if a
> Gauss-Laguerre will work well.
>
> K =@(t) exp(b0+b1.*t)./(1+exp(b0+b1.*t))./m;
> ezplot(K)
>
> It is not really that polynomial, although it is well behaved,
> with a sigmoidal shape.
>
> I then tried a Gauss-Laguerre quadrature on this function. It
> is not polynomial enough that a low order Gauss-Laguerre
> works, and roots fails above 14 nodes to generate any higher
> order Gauss-Laguerre rules than that. (I should have used a
> better routine to solve for the Gauss rule. No matter, since
> there are better choices to be found.) Look instead at a
> numerical tool like quadgk.
>
> quadgk(fun,0,30,'reltol',1e-12)
> ans =
> 0.584687911313577
>
> quadgk(fun,0,100,'reltol',1e-12)
> ans =
> 0.584733311243336
>
> quadgk(fun,0,200,'reltol',1e-12)
> ans =
> 0.584733311243339
>
> Or, try quadl. Note that quadl used more than 1200
> function evaluations to get this result, but I had set
> the tolerance rather small.
>
> [I,fev] = quadl(fun,0,200,1e-12)
> I =
> 0.584733311243339
> fev =
> 1248
>
> The default tolerance for quadl gave me a reasonable
> result at a cost of only 168 function evaluations. In fact
> this is far better accuracy than I would get from trapz
> anyway.
>
> [I,fev] = quadl(fun,0,200)
> I =
> 0.584733310572222
> fev =
> 168
>
> How well does trapz do? At a cost of 3000 function
> evaluations, you only get about 4 digits of accuracy,
> whereas the tools like quadl or quadgk will yield
> twice as many correct digits for vastly less work.
>
> t = 0:0.01:30;
> trapz(t,fun(t))
> ans =
> 0.584687862068707
>
> So I'm not at all sure why you have chosen to use
> trapz here, perhaps only that you did not know these
> other tools exist in MATLAB. Trapezoidal integration
> is rarely a good choice.
>
> Can we do better yet? Certainly so, if I knew that your
> function would always have the same essential character
> of a sigmoid multiplied by a negative exponential.
>
> I'll stop here for now, unless you wanted me to explore
> this last avenue more deeply. Ask again if so.
>
> John

Thanx again John,

It's me again.

You were right that I didn't know these other tools in MATLAB. I tried them though and with trial and error I think that quadgk is the right choise. Although it didn't work for the next:

setb0=-4.5313
setb1=2.7481
mhat=2.6470
pc=quadgk(@(gr) exp(setb0+setb1.*gr)./(1+exp(setb0+setb1.*gr)).*1/mhat.*exp(-gr./mhat),0,Inf,'reltol',1e-12)

For small changes of the values setb0 setb1 mhat it might work. As you can suppose setb0, setb1 mhat are random variables and I get a value of them in each iteration of the code. Then I try to integrate the above function which is the (sigmoidal*exponential distribution) .

To tell you the whole problem I am interested in solving the following integrals for whatever values of setb0 setb1, and mhat>0.

Integrals for the exponential distribution:

pc= quadgk(@(gr) exp(setb0+setb1*gr)/(1+exp(setb0+setb1*gr))*1/mhat*exp(-gr/mhat),0,Inf);

muc1= quadgk(@(gr) exp(setb0+setb1*gr)/((1+exp(setb0+setb1*gr))^2)*1/mhat*exp(-gr/mhat),0,Inf);

muc2= quadgk(@(gr) gr*exp(setb0+setb1*gr)/((1+exp(setb0+setb1*gr))^2)*1/mhat*exp(-gr/mhat),0,Inf);

As you can see, in all the above there is a product (something * exponential distribution).

The integrals for another distribution ( specifically the extended generalized gamma with p.d.f.=(1/gamma(g)).*(g.^g)*(l.^(a.*g)).*(x.^(a.*g-1)).*exp(-g.*((l.*x).^a)) ) where -Inf<a<Inf , l>0, g>0, x>0.

are the corresponding ones that come if one substitute the pdf of exponential with the pdf of the GG.

Any help for the above will be more than welcome...

Maybe if you can guide me through the example I gave you on top I might find me way through to the rest...

Thanx again for your time John!

Re: trapz, numerical integration

leo nidas — Sun, 07 Jun 2009 18:20:16 -0400

Well after some searching I think this is it:

A new quadrature routine for improper and oscillatory integrals
Applied Mathematics and Computation 189 (2007) 452–461

E. Sermutlu , H.T. Eyyuboglu

Re: trapz, numerical integration

Matt Fig — Sun, 07 Jun 2009 19:03:01 -0400

Why not use regular old Gauss legendre:

setb0=-4.5313
setb1=2.7481
mhat=2.6470
pc=gaussquad(@(gr) exp(setb0+setb1.*gr)./(1+exp(setb0+setb1.*gr)).*1/mhat.*exp(-gr./mhat),0,200)
pc =

0.551709492450232

b0=-5;
b1=3;
m=3
pc= quadgr(@(t) exp(b0+b1.*t)./(1+exp(b0+b1.*t)).*1/m.*exp(-t./m),0,200)

pc =
0.584733311243339

http://www.mathworks.com/matlabcentral/fileexchange/8679

Re: trapz, numerical integration

John D'Errico — Sun, 07 Jun 2009 21:44:02 -0400

"Matt Fig" <spamanon@yahoo.com> wrote in message <h0h2t5$juo$1@fred.mathworks.com>...
> Why not use regular old Gauss legendre:
>
> setb0=-4.5313
> setb1=2.7481
> mhat=2.6470
> pc=gaussquad(@(gr) exp(setb0+setb1.*gr)./(1+exp(setb0+setb1.*gr)).*1/mhat.*exp(-gr./mhat),0,200)
> pc =
>
> 0.551709492450232
>
>
>
> b0=-5;
> b1=3;
> m=3
> pc= quadgr(@(t) exp(b0+b1.*t)./(1+exp(b0+b1.*t)).*1/m.*exp(-t./m),0,200)
>
> pc =
> 0.584733311243339
>
>
>
> http://www.mathworks.com/matlabcentral/fileexchange/8679

The problem is that that you need to decide how far to
integrate out. Yes, 200 works here. It is not too large that
the integration fails, seeing only zeros, returning zero.
Too small of a value for the upper limit and it fails again,
missing too much of the area. So if you choose a reasonable
value for that upper limit, any code of reasonable quality
will work nicely - quadgk, quadgr, quadl, etc.

I'll suggest that a reasonable compromise is to pick some
value cheaply, such that the integration need not do too
much work.

b0=-5;
b1=3;
m=3

Here is the function of interest:

  exp(b0+b1.*t)./(1+exp(b0+b1.*t)).*exp(-t./m)/m

Do a transformation,

  s = t/m

and play around a bit, so the function becomes

  exp(-s)/(1 + exp(-(b0 + b1*m*s)))

or, this

  1/(exp(s) + exp(-b0 + (1 - b1*m)*s))

At s = 0, we get 1/(1+exp(-b0)). It will grow to some
maximum value, and then head back down to zero.

Where does the max occur? When the denominator is
minimized. This happens at the point:

   smax = (-b0 + log(b1*m - 1))/(b1*m)
smax =
         0.786604615742204

for these values of b0, b1, m.

quadgr(funs,0,100*smax,1.e-13)
ans =
         0.584733311243334

quadgk(funs,0,100*smax,'reltol',1.e-13)
ans =
         0.584733311243339

We can integrate out to some reasonably large multiple
of smax, and be confident that our result is correct.
Perhaps even better is to use a composite rule. We
know this function behaves like exp(-s) for large
enough s. 10*smax will surely be reasonably large.

So faster would seem to be to use a Gauss-Laguerre
rule beyond that point, adding the result to that of
your favorite code for small s. (In fact, quadgr seems
to be a bit faster here than quadgk.) The composite
method would seem to maximize both speed and
accuracy.

% ====================================
function [I,smax] = lapint(b0,b1,m)
% smax is the maximum of the integrand.
% the break point between quadgr and the
% laguerre will be 10*smax
smax = (-b0 + log(b1*m - 1))/(b1*m);

funs = @(s) 1./(exp(s) + exp(-b0 + (1-b1*m)*s));
I = quadgr(funs,0,10*smax,1e-12);

% now, use gauss-Laguerre for the rest, out to inf.
s0 = 10*smax;
Lagfun = @(u) 1./(1+exp(-b0-b1*m*(u+s0)));

nodes = [0.229680505425136, 0.772144910375415, 1.63105309906745, ...
          2.81514459001226, 4.33716407733756, 6.21464276455924, ...
          8.47116398134672, 11.1383319657508, 14.2589100216242, ...
          17.8920534381695, 22.1226201748336, 27.0793114990475, ...
          32.9749735523974, 40.2165837114966, 49.8462217085566];
weights = [0.0705024286637906, 0.249558909040495, 0.325533115492132, ...
         0.227727088139343, 0.0961880579830968, 0.0256342848941844, ...
         0.00435664720289306, 0.000467455769950909, 3.07992383129659e-05, ...
      1.18839046780406e-06, 2.49313557094408e-08, 2.52956064020393e-10, ...
      1.01977766969811e-12, 1.12207052360831e-15, 1.30933713763764e-19];

ILag = exp(-s0)*dot(Lagfun(nodes),weights);
I = I + ILag;
% ====================================

Is it reasonably accurate? Yes. As accurate as are the
others.

[I,smax] = lapint(b0,b1,m)
I =
         0.584733311243339
smax =
         0.786604615742204

How fast is it?

timeit(@() lapint(b0,b1,m))
ans =
       0.00337505370568894

timeit(@() quadgr(funs,0,100*smax))
ans =
       0.00418888466636987

timeit(@() quadgk(funs,0,100*smax))
ans =
       0.00491137823164936

HTH,
John