Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2004
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2004

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: UnitStep function leads to very difficult Integration

  • To: mathgroup at smc.vnet.net
  • Subject: [mg48120] Re: UnitStep function leads to very difficult Integration
  • From: astanoff_otez_ceci at yahoo.fr (astanoff)
  • Date: Fri, 14 May 2004 00:12:19 -0400 (EDT)
  • References: <c7ut4m$q9l$1@smc.vnet.net>
  • Sender: owner-wri-mathgroup at wolfram.com

Nathan Moore wrote:

> Hello all,

> I'm trying to evaluate expectation values from joint probability 
> distributions which I've had to define in a piecewise manner, ie

> 	norm = Integrate[G[1,x]G[2,x],{x,0,10}];
> 	<x^2> = Integrate[G[1,x]G[2,x] x^2,{x,0,10}] * norm

> where the G[k,x] is a set of polynomials defined piecewise on 
> intervals, {(0,2),(2,4),(4,6),etc}.

> I assume the most Mathematica-friendly way to define these functions is 
> with the UnitStep[] function.  This works to varying degree.  I was 
> able to check the probability normalization this way (integrating only 
> one G[k,x]) with an exact result.  Unfortunately though Mathematica 
> seems unable to find an exact experssion for the expectation when the 
> G[k,x] polynomials take sufficiently complex form.

> The integration output starts to look like,
> 	Integrate[ d^4 ( (d-3)UnitStep[3-d] UnitStep[d-1]/d + UnitStep[1 - d] 
> UnitStep[d])^2, {d,0,3}]
[...]
This is the way you could get exact rational values :

In[1]:=
stepIntegrate[f_(*expression*),
      d_(*variable*),
      cv_List(*critical values*),
      pd_(*polynomial degree*),
      incr_(*increment*),
      eps_(* epsilon to avoid boundary integration troubles*)
      ] := Module[{termIntegrate, th, ti},
      termIntegrate[{a_, b_}] :=
        With[{},
          t=Table[{x, f /. d -> x},{x, a+eps, b-eps, incr}];
          tfit = Fit[t,Table[d^n,{n,0,pd}],d] // Chop // Rationalize;
          {a,b,tfit}];
      th = Thread[{Drop[cv,-1], Rest[cv]}];
      ti = termIntegrate /@ th;
      Plus @@ (Integrate[#[[3]],{d, #[[1]]+eps, #[[2]]-eps}]& /@ ti)
      ];

(* your example : *)
expr=d^4 ((d-3)UnitStep[3-d] UnitStep[d-1]/d+UnitStep[1-d] UnitStep[d])^2;
stepIntegrate[expr,d,{0,1,3},5,0.1,10.^-15]//Rationalize
Out[1]=33/5
--
0% de pub! Que du bonheur et des vrais adhérents !
Vous aussi inscrivez-vous sans plus tarder!!
Message posté à partir de http://www.gyptis.org, BBS actif depuis 1995.




  • Prev by Date: nonlinear theory of elasticity question using Mathematica
  • Next by Date: Maen
  • Previous by thread: UnitStep function leads to very difficult Integration
  • Next by thread: forward link : Re: Mandelbrot Set & Mathematica