Help Improving this integral calculation / solution

Help Improving this integral calculation / solution

Post by DOD » Mon, 08 Mar 2010 18:04:52


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?
 
 
 

Help Improving this integral calculation / solution

Post by Bob Hanlo » Tue, 09 Mar 2010 20:13:00


d = 7/10;
v = 5/100;
a = 1/10;
nk = (a + (1 - a) k v t)/(d + a + (1 - a) k v t) //
Simplify

(9*k*t + 20)/(9*k*t + 160)

pg[k_] = Exp[1 - k];

gavg = Integrate[k pg[k], {k, 1, Infinity}]

2

result = 1/gavg Integrate[nk k pg[k],
{k, 1, Infinity}, Assumptions -> 0 < t < 1]

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

test = FullSimplify[result]

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

root = FindRoot[test == t, {t, .1}]

{t -> 0.1420452943881465}

Plot[{t, test}, {t, 0, .3},
Epilog -> {Red, AbsolutePointSize[4],
Point[{t, t} /. root]}]


test - t /. root

-4.973799150320701*^-14


Bob Hanlon



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

 
 
 

Help Improving this integral calculation / solution

Post by danl » Tue, 09 Mar 2010 20:15:05

> I have an integral I need to calculate, with one variable left

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
 
 
 

Help Improving this integral calculation / solution

Post by David Par » Tue, 09 Mar 2010 20:15:37

This goes much better if you use exact values:

d = 7/10;
v = 1/20;
a = 1/10;
nk = (a + (1 - a) k v t)/(d + a + (1 - a) k v t) // Simplify

(20 + 9 k t)/(160 + 9 k 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];
f[t_] = result

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

Plot[f[t] - t, {t, 0, 0.2},
PlotRange -> Automatic]

FindRoot[f[t] == t, {t, .15}]
{t -> 0.142045}

I tried solving the integral with a, d and v as parameters, and accepting
the conditions but that did not work.


David Park
XXXX@XXXXX.COM
http://www.yqcomputer.com/ ~djmpark/


From: DOD [mailto: XXXX@XXXXX.COM ]


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?
 
 
 

Help Improving this integral calculation / solution

Post by DrMajorBo » Tue, 09 Mar 2010 20:16:10

Try this:

{d, v, a} = Rationalize@{.7, .05, .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]}]

2

result = 1/gavg Integrate[nk k pg[k], {k, 1, \[Infinity]},
Assumptions -> 0 < t < 1]

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

test = FullSimplify[result]

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

FindRoot[test == t, {t, .1}]

{t -> 0.142045}

Bobby





--
XXXX@XXXXX.COM
 
 
 

Help Improving this integral calculation / solution

Post by dh » Tue, 09 Mar 2010 20:41:41

Hi,
the variable "test" contains differences of nearly identical huges
terms. Therefore, machine accuracy is not good enough, your results are
simply calculation errors. You must increase the accuracy.
Daniel





--

Daniel Huber
Metrohm Ltd.
Oberdorfstr. 68
CH-9100 Herisau
Tel. +41 71 353 8585, Fax +41 71 353 8907
E-Mail:<mailto: XXXX@XXXXX.COM >
Internet:< http://www.yqcomputer.com/ >
 
 
 

Help Improving this integral calculation / solution

Post by DOD » Wed, 10 Mar 2010 20:20:18

Thanks for all the tips and explanations. As shown, these methods
work. I actually finally got decent results by doing the integration/
root solving simultaneously:

---
result=NMinimize[{ (1/gavg Integrate[nk k pg[k], {k, 1, \[Infinity]},
Assumptions -> 0 < t < 1] - t)^2,0<t<1} , t];
answer = t/.result[[2]]
---

And this found the t I wanted. This was after haphazard reformulation
of the problem into various forms amenable to various built-in
functions Thanks for these methods and insight into what the problem
was; I need to do this calculation for various parameters (d,v,a), so
any increase in speed is key.