MathGroup Archive 2010

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

Search the Archive

Re: Help Improving this integral calculation / solution

  • To: mathgroup at smc.vnet.net
  • Subject: [mg108100] Re: [mg108044] Help Improving this integral calculation / solution
  • From: danl at wolfram.com
  • Date: Mon, 8 Mar 2010 06:14:44 -0500 (EST)
  • References: <201003070903.EAA00857@smc.vnet.net>

> I have an integral I need to calculate, with one variable left
> symbolic, and then use the result to find a numerical solution to an
> equation, and the result of the integration is not giving me what I
> need:
> ------
> d = .7;
> v = .05;
> a = .1;
> nk = (a + (1 - a) k v t)/(d + a + (1 - a) k v t);
>
> pg[k_] = Exp[1-k];
> gavg = Integrate[k pg[k], {k, 1, \[Infinity]}]; (* This is just 2 *)
> result =   1/gavg Integrate[nk k pg[k], {k, 1, \[Infinity]},
> Assumptions -> 0 < t < 1];
> test = FullSimplify[result]
>
> FindRoot[test == t, {t, .1}]
>
> -----
> This code always gives up and stays at the initial guess.  So, I look
> at the result of the integration, "result" (or it's simplified
> version, test2) and calculated it for various values of t, and it is
> always zero.
> ----
> test/.t->{.1,.2,.3}
> ----
> Output:{3.36999*10^66, 0., 0.}
> -----
>
> So that's a problem.  If I set t=.3, say, and the beginning, and
> calculate the integral, I get 0.160047, which is clearly not zero.
> -----
> d = .7;
> v = .05;
> a = .1;
> t = .3;
> nk = (a + (1 - a) k v t)/(d + a + (1 - a) k v t);
>
> pg[k_] = Exp[1 - k];
> gavg =Integrate[k pg[k], {k, 1, \[Infinity]}];(*This is just 2*)
>
> result = 1/gavg Integrate[nk k pg[k], {k, 1, \[Infinity]}]
> ----
> Output=0.160047
> -----
>
>  So there is a problem in the Integrate step.  So I want to find a
> point where the output of that integral (result), as a function of t,
> satisfies result=t.  Is there anyway to do this using built-in
> functions?

One simple method is to use NIntegrate, restricted so it does not try to
evaluate for nonnumeric input.

res[t_?NumberQ] := NIntegrate[nk k pg[k], {k, 1, \[Infinity]}]/gavg

In[150]:= FindRoot[res[t] == t, {t, .1}]
Out[150]= {t -> 0.142045}

Can also do this using the original approach with Integrate. In that case
I recommend using exact values for the parameters.

In[159]:= d = 7/10;
v = 1/20;
a = 1/10;
nk = (a + (1 - a) k v t)/(d + a + (1 - a) k v t);
pg[k_] = Exp[1 - k];
gavg = Integrate[k pg[k], {k, 1, \[Infinity]}];
result = 1/gavg Integrate[nk k pg[k], {k, 1, \[Infinity]},
   Assumptions -> 0 < t < 1]
FindRoot[result == t, {t, .1}]

Out[165]= (
9 t (-70 + 9 t) -
 11200 E^(1 + 160/(9 t)) ExpIntegralEi[-1 - 160/(9 t)])/(81 t^2)

Out[166]= {t -> 0.142045}

Daniel Lichtblau
Wolfram Research




  • Prev by Date: Re: Possible bug in Table and/or lists (Mathematica v.7)
  • Next by Date: Re: Re: Write an expression in a specific form
  • Previous by thread: Re: Help Improving this integral calculation / solution
  • Next by thread: Re: Help Improving this integral calculation / solution