Solving for upper bound of integral with known value

Question

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

8. +(-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]

Answer 1

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;

Compare with @Bob's answer:

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

Answer 2

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 *)