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 21:44:02 +0000 (UTC)
Organization: John D'Errico (1-3LEW5R)
Lines: 146
Message-ID: <h0hcb2$ant$1@fred.mathworks.com>
References: <h048vl$2v$1@fred.mathworks.com> <h04bpo$2vo$1@fred.mathworks.com> <h0fudh$pfo$1@fred.mathworks.com> <h0gade$kfk$1@fred.mathworks.com> <h0gr9d$6db$1@fred.mathworks.com> <h0h0d0$897$1@fred.mathworks.com> <h0h2t5$juo$1@fred.mathworks.com>
Reply-To: "John D'Errico" <woodchips@rochester.rr.com>
NNTP-Posting-Host: webapp-05-blr.mathworks.com
Content-Type: text/plain; charset="ISO-8859-1"
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1244411042 11005 172.30.248.35 (7 Jun 2009 21:44:02 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Sun, 7 Jun 2009 21:44:02 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 869215
Xref: news.mathworks.com comp.soft-sys.matlab:545382


"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