Path: news.mathworks.com!not-for-mail
From: "John D'Errico" <woodchips@rochester.rr.com>
Newsgroups: comp.soft-sys.matlab
Subject: Re: trapz, numerical integration
Date: Sun, 7 Jun 2009 12:05:02 +0000 (UTC)
Organization: John D'Errico (1-3LEW5R)
Lines: 110
Message-ID: <h0gade$kfk$1@fred.mathworks.com>
References: <h048vl$2v$1@fred.mathworks.com> <h04bpo$2vo$1@fred.mathworks.com> <h0fudh$pfo$1@fred.mathworks.com>
Reply-To: "John D'Errico" <woodchips@rochester.rr.com>
NNTP-Posting-Host: webapp-02-blr.mathworks.com
Content-Type: text/plain; charset="ISO-8859-1"
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1244376302 20980 172.30.248.37 (7 Jun 2009 12:05:02 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Sun, 7 Jun 2009 12:05:02 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 869215
Xref: news.mathworks.com comp.soft-sys.matlab:545316


"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