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.