Problems with a constant term in an Integrate command

Question

When asking Mathematica to integrate a function with Integrate, I'm getting drastically different results depending on whether I have a proportionality constant inside or outside of the Integrate command.

To summarize: I'm attempting to define a probability distribution from the function ds:

ds[e_, th_]:= a^2r^2(P[e, th]^2)(P[e, th]+P[e, th]^(-1)-1+Cos[th]^2)/2

with

P[e_, th_]:=1/(1+(e/(m*c^2))(1-Cos[th]))

and a, r, c, and m real constants, as follows:

c = 3*10^(8);
h = 1.05*10^(-34);
m = 9.11*10^(-31);
r = h/(m*c);
a = 7.30*10^(-3);

I'd like the distribution to give the probability density function for values of th from 0 to 2pi. However, for this to work properly with ProbabilityDistribution, ds needs to be normalized. I should be able to do this by integrating ds over the relevant interval (for a specific e, in this case 10^(-15)) and dividing ds by the resulting constant:

s=Integrate[ds[10^(-15), th],{th, 0, 2 Pi}]

1.66118*10^(-25)

ds2[th_]:=ds[10^(-15), th]/s

However, something seems to be broken in this solution. If I define a probability distribution in terms of ds2, it doesn't act as though the distribution is normalized, and instead returns values far outside the desired range. To diagnose the problem, I've tried integrating ds divided by s. Here's what I'm getting: When s is outside the integral, I have

Integrate[ds[10^(-15), th], {th, 0, 2 Pi}]/s

1.

as expected, but when s is inside the integral, I get

Integrate[ds[10^(-15), th]/s, {th, 0, 2 Pi}]

0.000217878

Because s is a constant, moving it inside the integral shouldn't affect the result, so clearly something is going wrong. I've retried the computation with common functions instead of something complicated like ds, and don't have any problems there. What might the problem be?

Answer 1

You are having round-off problems. One way to eliminate them is by replacing your approximate real numbers by rational numbers, because Mathematica computes results with exact numbers to arbitrary precison.

c = 3*10^(8); h = 105*10^(-36); m = 911*10^(-33); r = h/(m*c); a = 730*10^(-5);

Then, both normalized integrals yield precisely 1.