[Top][All Lists]

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

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

From: Vasu Jaganath
Subject: Re: [Help-gsl] fixed point or adaptive integration for calculating moments using beta PDF?
Date: Sun, 31 Dec 2017 17:28:36 -0700

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

Attachment: test_snippet
Description: Binary data

reply via email to

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