# Solving for upper bound of integral with known value

I'm trying to calculate the upper bound of an integral. My current code is

f[r_?NumericQ] =
(1/(Pi*a^3))*4*Pi*(r^2)*E^(-r/a);
Plot[f[r], {r, 0, 1}]


When I evaluate this, the graph is incorrect. The current graph and error I'm getting are

last night, however, I was able to get the graph to look like this

What changed? Why isn't the current code working? Alternatively, are there other ways to do this calculation?

EDIT: Using the code given I redid my calculations. The new code is

a = 5.291*10^-11;

f[r_] = (1/(Pi*a^3))*4*Pi*(r^2)*E^(-r/a);

Plot[f[r], {r, 0, 1}]
clear[b]
g[b_] = Integrate[f[r], {r, 0, b}]
sol = g[b] /. sol
NSolve[{g[b] == 1/2, b > 0}, b][[1]]


the output is

1. +(-8.-1.512*10^11 b-1.42884*10^21 b^2) E^(-1.89*10^10 b) ReplaceAll::reps: {sol} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing. \$RecursionLimit::reclim2: Recursion depth of 1024 exceeded during evaluation of >8. +(-8.-1.512*10^11 b-1.42884*10^21 b^2) E^(-1.89*10^10 b). Hold[8. +(-8.-1.512*10^11 b-1.42884*10^21 b^2) E^(-1.89*10^10 b)/. sol]
• you need a and b constants – J42161217 Oct 14 '18 at 16:30
• Please use an informative title, that will help others find solutions to similar problems. – David G. Stork Oct 14 '18 at 16:42

Your problem can be solved exactly. The integral is:

f[r_] := (1/(Pi*a^3))*4*Pi*(r^2)*E^(-r/a)

int = Integrate[f[r], {r, 0, b}]


8 - (4 (2 a^2 + 2 a b + b^2) E^(-(b/a)))/a^2

Let $$b = s a$$:

scaledint = int /. b->a s //Simplify


8 - 4 E^-s (2 + 2 s + s^2)

Set this result to $$1/2$$ and solve for s:

rt = s /. First @ Solve[scaledint == 1/2 && s>0, s]


Root[{16 - 15 E^#1 + 16 #1 + 8 #1^2 &, 0.89782942356979330033}]

So, the dependence of b on a is:

b[a_] = a rt;


b[1/2]
%//N


1/2 Root[{16 - 15 E^#1 + 16 #1 + 8 #1^2 &, 0.89782942356979330033}]

0.448915

and:

b[5291 10^-14]
%//N


(5291 Root[{16 - 15 E^#1 + 16 #1 + 8 #1^2 &, 0.89782942356979330033}])/100000000000000

4.75042*10^-11

Since the definition of f does not use numeric techniques (e.g., NSolve, FindRoot, NIntegrate) , there is no reason to restrict its argument to being numeric by using NumericQ. However, in order to Plot, the parameter a must have a numeric value which you presumably defined but did not provide.

It appears that you used

a = 1/2;

f[r_] = (1/(Pi*a^3))*4*Pi*(r^2)*E^(-r/a);

Plot[f[r], {r, 0, 1}]


Clear[b]


f can be integrated analytically so use of NIntegrate is unnecessary.

g[b_] = Integrate[f[r], {r, 0, b}]

(* 8 + 8 (-1 - 2 b (1 + b)) E^(-2 b) *)


The equation g[b] == 1/2 can be solved using NSolve provided that you restrict b to be positive. This avoids use of FindRoot and the need to provide an initial estimate.

sol = NSolve[{g[b] == 1/2, b > 0}, b][[1]]

(* {b -> 0.448915} *)


Verifying

g[b] /. sol

(* 0.5 *)


EDIT: For the revised value of a

a = 5291*10^-14;

f[r_] = (1/(Pi*a^3))*4*Pi*(r^2)*E^(-r/a);


For a with such a small value, the values of E^(-r/a) are too small to evaluate. Plot the Log of f

logf[r_] = Log[f[r]] // PowerExpand // FullSimplify

(* -((100000000000000 r)/5291) + 44 Log[2] + 42 Log[5] - 3 Log[5291] + 2 Log[r] *)

Plot[logf[r], {r, 0, 1}]


Clear[b]
g[b_] = Integrate[f[r], {r, 0, b}]

(* 8 - (1/27994681)
8 (27994681 +
100000000000000 b (5291 + 50000000000000 b)) E^(-100000000000000 b/5291) *)


You did not copy the code correctly. If done properly,

sol = NSolve[{g[b] == 1/2, b > 0}, b][[1]]

(* {b -> 4.75042*10^-11} *)


Verifying

g[b] /. sol

(* 0.5 *)

• When I put in your suggestions, it removed the error message but now the solutions I get include b in them as opposed to giving me an output of a value for b. – Dom Oct 14 '18 at 17:27
• If you use the code that I posted you should get the results that I show. If not, start with a fresh kernel. If after a restart you still are having problems post the code (not an image). – Bob Hanlon Oct 14 '18 at 17:33