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