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: Vasu Jaganath
Subject: Re: [Help-gsl] fixed point or adaptive integration for calculating moments using beta PDF?
Date: Mon, 01 Jan 2018 04:12:30 +0000

That's what I am trying to do, the tables are in a weird format

On Sun, Dec 31, 2017, 9:09 PM Patrick Alken <address@hidden> wrote:

> 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> 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
>>
>> 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
>> >>>
>> >>
>>
>>
>>
>
>


reply via email to

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