help-gsl
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [Help-gsl] fixed point or adaptive integration for calculating momen


From: Patrick Alken
Subject: Re: [Help-gsl] fixed point or adaptive integration for calculating moments using beta PDF?
Date: Sun, 31 Dec 2017 21:09:10 -0700
User-agent: Mozilla/5.0 (X11; Linux x86_64; rv:52.0) Gecko/20100101 Thunderbird/52.5.0

Why not just plot Q(x) over your integration interval?

On 12/31/2017 09:00 PM, Vasu Jaganath wrote:
> As far as I know, my Q doesn't contain any singularities but I will
> recheck my tables again just to confirm.
>
> Thanks,
> Vasu
>
> On Sun, Dec 31, 2017 at 8:37 PM, Patrick Alken <address@hidden
> <mailto:address@hidden>> 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#fixed-point-quadratures
>     
> <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
>     <address@hidden <mailto:address@hidden>> wrote:
>     >
>     >> Can you give a concrete example of a typical function Q?
>     >>
>     >> On Sat, Dec 30, 2017 at 3:42 AM, Vasu Jaganath
>     <address@hidden <mailto:address@hidden>>
>     >> 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
>     >>>
>     >>
>
>
>



reply via email to

[Prev in Thread] Current Thread [Next in Thread]