0
$\begingroup$

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 error graph

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

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]
$\endgroup$
  • 1
    $\begingroup$ you need a and b constants $\endgroup$ – J42161217 Oct 14 '18 at 16:30
  • $\begingroup$ Please use an informative title, that will help others find solutions to similar problems. $\endgroup$ – David G. Stork Oct 14 '18 at 16:42
2
$\begingroup$

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

$\endgroup$
1
$\begingroup$

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}]

enter image description here

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}]

enter image description here

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 *)
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Dom Oct 14 '18 at 17:27
  • $\begingroup$ 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). $\endgroup$ – Bob Hanlon Oct 14 '18 at 17:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.