Changing constants is affecting integration

When I define a function in the following way, and try evaluating it at specific values of its arguments, Mathematica returns an answer.

xa[δx_, δy_, e1_, e2_, W_, L_, r_] := NIntegrate[
2*e1*e2*Exp[-((x - δx)^2 + (y - δy)^2)/W^2]*Exp[-(x^2 + y^2)/L^2]
, {x, 0, r}
, {y, -(r^2 - x^2)^(1/2), (r^2 - x^2)^(1/2)}
, MinRecursion -> 10
, MaxRecursion -> 20
, PrecisionGoal -> 3
, WorkingPrecision -> 40
]

However if I modify the function definition slightly by merely replacing the argument e2 by its expression in terms of another argument L (which are both constants with respect to this integration), I do not get any answer.

xa1[δx_, δy_, e1_, W_, L_, r_] := NIntegrate[
2*e1*(L^-1*Sqrt[P/(π*ϵ0*c)])*Exp[-((x - δx)^2 + (y - δy)^2)/W^2]*Exp[-(x^2 + y^2)/L^2]
, {x, 0, r}
, {y, -(r^2 - x^2)^(1/2), (r^2 - x^2)^(1/2)}
, MinRecursion -> 10
, MaxRecursion -> 20
, PrecisionGoal -> 3
, WorkingPrecision -> 40
]

Why does this happen? How could I make the second function work? Thanks.

• Can you provide an example of the arguments you are using to call these functions, both for the case it works and the case it doesn't? I have formatted your code for readability, please check that there are no differences with your code. – rhermans Sep 27 '17 at 15:06
• Is P defined? What is it? – John Doty Sep 27 '17 at 15:09

Let's factor the constants out of NIntegrate and see what progress we can make. In the following we re-define the two functions with constants multiplying the integral, we calculate the ratio of the two functions (ignoring error messages) and find the required value of P in terms of e2, and vice versa:

ClearAll["Global`*"]
xa[δx_, δy_, e1_, e2_, W_, L_, r_] :=

2*e1*e2*NIntegrate[
Exp[-((x - δx)^2 + (y - δy)^2)/W^2]*
Exp[-(x^2 + y^2)/L^2], {x, 0,
r}, {y, -(r^2 - x^2)^(1/2), (r^2 - x^2)^(1/2)},
MinRecursion -> 10, MaxRecursion -> 20, PrecisionGoal -> 3,
WorkingPrecision -> 40]

xa1[δx_, δy_, e1_, W_, L_, r_] :=
2*e1*(L^-1*Sqrt[P/(π*ϵ0*c)])*
NIntegrate[
Exp[-((x - δx)^2 + (y - δy)^2)/W^2]*
Exp[-(x^2 + y^2)/L^2], {x, 0,
r}, {y, -(r^2 - x^2)^(1/2), (r^2 - x^2)^(1/2)},
MinRecursion -> 10, MaxRecursion -> 20, PrecisionGoal -> 3,
WorkingPrecision -> 40]

ratio = xa[δx, δy, e1, e2, W, L, r]/
xa1[δx, δy, e1, W, L, r]
Solve[ratio == 1, P] // Flatten
Solve[ratio == 1, e2] // Flatten

(*  (e2 L Sqrt[π])/Sqrt[P/(c ϵ0)]

{P -> c e2^2 L^2 π ϵ0}

{e2 -> Sqrt[P/(c ϵ0)]/(L Sqrt[π])}  *)

Note that we are passing symbols to our functions so we must ignore the error messages that NIntegrate gives. Since we are calculating the ratio of the two functions and the only difference in the two functions is the constants that are outside the integrals, Mathematica can reduce the ratio as shown.

It is true that we should get the same thing with the constants inside the NIntegrate, but it is easier to debug with them outside.

But to answer your question, we can make the second function work by setting P = c e2^2 L^2 π ϵ0 before we define it.

Also, as a check, we use our calculated value of e2 as follows, again ignoring the NIntegrate warnings:

xa[δx, δy, e1, Sqrt[P/(c ϵ0)]/(L Sqrt[π]), W, L, r] ==
xa1[δx, δy, e1, W, L, r]

(*  True  *)