[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 20:37:23 -0700 |
User-agent: |
Mozilla/5.0 (X11; Linux x86_64; rv:52.0) Gecko/20100101 Thunderbird/52.5.0 |
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
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> 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>
>> 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
>>>
>>