You might be able to used fixed quadrature - you will probably need a large number of nodes to capture that sharp feature. I would compare both the CQUAD and fixed-point methods to see which works better for you.
Patrick On 12/31/2017 11:10 PM, Vasu Jaganath wrote: > HI Martin, Yes one of my Q is very discontinuous with respect to my > integration variable Z. > > I have attached the plots for you to see > > On Sun, Dec 31, 2017 at 9:20 PM, Vasu Jaganath > <[email protected] <mailto:[email protected]>> wrote: > > Yes, I will show you the plots soon, Q is actually 2 variable > function but for my calculations I am treating one of the > variables as a parameter, which is a physically valid assumption. > Yes I do encounter alpha and beta values less than 1. > > > On Sun, Dec 31, 2017, 9:13 PM Martin Jansche <[email protected] > <mailto:[email protected]>> wrote: > > So you want to find E[f] = \int_0^1 f(x) dbeta(x | a, b) dx. > Can you plot > your typical f? And what are typical values of the parameters > a and b? Do > you encounter a<=1 or b<=1? If so, how does f(x) behave as x > approaches 0 > or 1? > > On Mon, Jan 1, 2018 at 3:37 AM, Patrick Alken > <[email protected] <mailto:[email protected]>> wrote: > > > The question is whether your Q contains any singularities, > or is highly > > oscillatory? Is such cases fixed point quadrature generally > doesn't do > > well. If Q varies fairly smoothly over your interval, you > should give > > fixed point quadrature a try and report back if it works > well enough for > > your problem. The routines you want are documented here: > > > > http://www.gnu.org/software/gsl/doc/html/integration.html# > <http://www.gnu.org/software/gsl/doc/html/integration.html#> > > fixed-point-quadratures > > > > Also, if QAGS isn't working well for you, try also the CQUAD > routines. > > I've found CQUAD is more robust than QAGS in some cases > > > > On 12/31/2017 05:28 PM, Vasu Jaganath wrote: > > > I have attached my entire betaIntegrand function. It is a > bit complicated > > > and very boiler-plate, It's OpenFOAM code (where scalar = > double etc.) I > > > hope you can get the jist from it. > > > I can integrate the Q using monte-carlo sampling integration. > > > > > > Q is nothing but tabulated values of p,rho, mu etc. I > lookup Q using the > > > object "solver" in the snippet. > > > > > > I have verified evaluating <Q> when I am not using those > <Q> values back > > in > > > the solution, It works OK. > > > > > > Please ask me anything if it seems unclear. > > > > > > > > > > > > > > > > > > > > > On Sun, Dec 31, 2017 at 3:55 PM, Martin Jansche > <[email protected] <mailto:[email protected]>> > > wrote: > > > > > >> Can you give a concrete example of a typical function Q? > > >> > > >> On Sat, Dec 30, 2017 at 3:42 AM, Vasu Jaganath < > > [email protected] <mailto:[email protected]>> > > >> wrote: > > >> > > >>> Hi forum, > > >>> > > >>> I am trying to integrate moments, basically first order > moments <Q>, > > i.e. > > >>> averages of some flow fields like temperature, density > and mu. I am > > >>> assuming they distributed according to beta-PDF which is > parameterized > > on > > >>> variable Z, whose mean and variance i am calculating > separately and > > using > > >>> it to define the shape of the beta-PDF, I have a cut off > of not using > > the > > >>> beta-PDF when my mean Z value, i.e <Z> is less than a > threshold. > > >>> > > >>> I am using qags, the adaptive integration routine to > calculate the > > moment > > >>> integral, however I am restricted to threshold of <Z> = > 1e-2. > > >>> > > >>> It complains that the integral is too slowly convergent. > However > > >>> physically > > >>> my threshold should be around 5e-5 atleast. > > >>> > > >>> I can integrate these moments with threshold upto 5e-6, > using > > Monte-Carlo > > >>> integration, by generating random numbers which are > beta-distributed. > > >>> > > >>> Should I look into fixed point integration routines? > What routines > > would > > >>> you suggest? > > >>> > > >>> Here is the (very simplified) code snippet where, I > calculate alpha and > > >>> beta parameter of the PDF > > >>> > > >>> zvar = min(zvar,0.9999*zvar_lim); > > >>> alpha = zmean*((zmean*(1 - > zmean)/zvar) - 1); > > >>> beta = (1 - zmean)*alpha/zmean; > > >>> > > >>> // inside the fucntion to be integrated > > >>> // lots of boilerplate for Q(x) > > >>> f = Q(x) * gsl_ran_beta_pdf(x, > alpha, beta); > > >>> > > >>> // my integration call > > >>> > > >>> helper::gsl_integration_qags (&F, 0, > 1, 0, 1e-2, > > 1000, > > >>> w, > &result, &error); > > >>> > > >>> And also, I had to give relative error pretty large, > 1e-2. However > > some of > > >>> Qs are pretty big order of 1e6. > > >>> > > >>> Thanks, > > >>> Vasu > > >>> > > >> > > > > > > > >
