[?] Use the solution of FindRoot as a integrate limit?

Questions
, , ,
1

For example, I want to do a numerical integrate as follow:

NIntegrate[1, {gamma, 1, 2}, {t, 10^6, 
  t /. FindRoot[
     1/21 t^(19/116) (116 B (t^(21/116) - t0^(21/116))) == 1/
      gamma, {t, 10^6}][[1]]}]

Here B and t0 are constants.

My goal is I want to do a multiple integrate, but the upper limit of the first integrate is the function of the second variable for example, the upper limit t(gamma) is the solution of equation:

1/21 t^(19/116) (116 B (t^(21/116) - t0^(21/116))) == 1/ gamma

so firstly I need to solve the equation numerically for gamma and find the value of upper limit, then do the integrate.

But the result returns:

FindRoot::nlnum: The function value {0. -1./gamma} is not a list of numbers with dimensions {1} at {t} = {1.*10^6}.
ReplaceAll::reps: {0.00608348 (-10 10^(5/58)+t^(21/116)) t^(19/116)==1/gamma} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing.
ReplaceAll::reps: {0.00608348 (-10 10^(5/58)+t^(21/116)) t^(19/116)==1/gamma} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing.
ReplaceAll::reps: {0.00608348 (-12.1957+t^(21/116)) t^(19/116)==1/gamma} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing.
General::stop: Further output of ReplaceAll::reps will be suppressed during this calculation.
NIntegrate::nlim: t = t/. 0.00608348 (-10 10^(5/58)+t^(21/116)) t^(19/116)==1/gamma is not a valid limit of integration.

Does anyone know what's wrong in my cod

e? And how to achieve my goal? Thank you!

2

In a first step, you should always try to get an intuition about what you are trying to calculate. Mathematica can help you, but it cannot replace thinking about your problem. You didn’t provide values for your constants b and t0, but let me assume they are greater zero. In a first step, you could have created a small manipulate showing the function whose root your are trying to find

Manipulate[
 Plot[
  1/21 t^(19/116) (116 b (t^(21/116) - t0^(21/116))) - 1/gamma, {t, 0,
    10}],
 {gamma, 1, 2},
 {t0, 0, 10},
 {b, 0, 10}
 ]

Mathematica graphics

When you move the sliders, then you see that as long as both constants are greater zero, you will always have a root at the very beginning. This lets me question your choice of the initial value for FindRoot which you chose to be 10^6. I guess with a value of 10, the root can be found fast.

In a next step, you should always ensure that your root-finding indeed returns a numerical value. NIntegrate will try to evaluate your FindRoot analytically, which will end badly. Therefore, define a function that ensures it can only be called for numeric values and returns the root you want to use as limit in your integration

root[gamma_?NumericQ] := Block[{t}, 
  t /. FindRoot[
    1/21 t^(19/116) (116 b (t^(21/116) - t0^(21/116))) == 1/gamma, {t, 10}]
  ]

As you see in the plot above, your root will be smaller than 10^6. I’m not sure you really want to integrate from 10^6 to e.g. 2 or the other way around. Let me assume you want it the other way around. If you have no experience with integrals, the easy way, in this case, is to solve it analytically but you leave out your real root function and replace it with an arbitrary function f. Then, your result takes the following form

Integrate[1, {gamma, 1, 2}, {t, f[gamma], 10^6}]
\int_1^2 (1000000-f(\gamma )) \, d\gamma

With integration rules, this simplifies to

10^6 - \int_1^2 f(\gamma ) \, d\gamma

where the function f is your root function. Let’s write this down:

10^6 - NIntegrate[root[\[Gamma]], {\[Gamma], 1, 2}]

This gives 999997. If you like, you can of course use your original call now:

NIntegrate[1, {gamma, 1, 2}, {t, root[gamma], 10^6}]

which gives the same result. Remember that I turned root[gamma] and 10^6 around! I’m not sure I deciphered exactly what you want but the points I made should help you in any case.

3

Thank you for your explanation!

Here the constants are B=0.00110132, t0=10^6, then you will find the root is larger than 10^6. And let’s assume the function I want to integrate is much more complicated than 1 so that we don’t have a analytical solution.

If you need the original function, it’s

gamma^(-2.2)(t0/t)^(57/290)(1 -
gamma116/21Bt^(19/116)(t^(21/116) - t0^(21/116)))^(1/5)

So should I just try your method:

root[gamma_?NumericQ] := Block[{t},
t /. FindRoot[
1/21 t^(19/116) (116 b (t^(21/116) - t0^(21/116))) == 1/gamma, {t, 10}]
]

so I can get a numerical value returned, then I can do my integrate?

4

Here the constants are B=0.00110132, t0=10^6, then you will find the root is larger than 10^6.

Excellent, then your insight is correct and I indeed used the wrong order of arguments.

I want to integrate is much more complicated than 1 so that we don’t have an analytical solution

Yes, you can essentially use my exact approach. Since your roots are indeed greater then 10^6, you should set your initial point in root back to 10^6 and you are ready to go:

b = 0.00110132;
t0 = 10^6;

root[gamma_?NumericQ] := 
 Block[{t}, 
  t /. FindRoot[1/21 t^(19/116) (116 b (t^(21/116) - t0^(21/116))) == 1/gamma, {t, 10^6}]
]

NIntegrate[
 gamma^(-2.2) (t0/t)^(57/290) (1 -  gamma*116/21 b*t^(19/116) (t^(21/116) - t0^(21/116)))^(1/5),
 {gamma, 1, 2}, {t, 10^6, root[gamma]}
]

This gives 4.11988*10^6 on my machine.

Btw, please don’t use variables (especially one letter variables) that start with a capital letter. You would be surprised how many of them are already used by Mathematica itself and you are likely to run into problems soon:

Select[Names["System`*"], StringMatchQ[#, CharacterRange["A", "Z"]] &]

(* {"C", "D", "E", "I", "K", "N", "O"} *)

Btw #2: If you like to participate more often here, you should lookup how Markdown syntax. It is an easy way to mark code as code (so that it is highlighted as in my post) and includes many more features like bold, italic, headlines, links, images, etc..