I am working with a function of the form:
import sympy as sm
a,b,x = sm.symbols('a b x')
f = sm.integrate(sm.cos(a*x + 0.5*b*x**2),(x,0,1))
The full integral looks like this in SymPy's formatted output. I don't really know the mathematical theory behind it, but I assume that this is a correct formulation.
However, the problem I'm having occurs in the arguments of the Fresnel S and C functions, where there is a factor of 1/sqrt(b). This means that when b is negative, the denominator becomes complex so the total argument discontinuously switches sign. The integral itself is always real as the complex terms cancel, but the sign change propagates to the overall value of the integral as well, eg:
f.subs([(a,1),(b,1e-9)]).evalf()
# =0.84147
f.subs([(a,1),(b,-1e-9)]).evalf()
# =−0.84147
At least for my application, it would make sense that this function is continuous at b=0. I am also a bit confused around the ambiguity of this reciprocal square root term, which is seems could be equivalently written: 1/sqrt(b) = sqrt(1)/sqrt(b) = sqrt(1/b), which loses the sign change for negative b. In fact my current workaround is to rewrite the Fresnel arguments from 1.0*a/(sqrt(pi)*sqrt(b)) to 1.0*a/(sqrt(pi))*sqrt(1/b) in the SymPy expression.
Additionally, I have tried this in MATLAB's symbolic toolbox, and it seems to give the output I am expecting by default:
syms a b x real
f = int(cos(a*x + 0.5*b*x^2),x,0,1)
MATLAB's formatted integral looks like this
Notably, the 1/sqrt(b) term appears as sqrt(1/b) in the numerator of the S and C arguments (sigma1 and sigma4 replacements).
How can I get SymPy to evaluate this function the same way MATLAB does? As mentioned, I have a hacky workaround of manually modifying part of the SymPy string expression, but this is obviously not ideal and probably won't work for more complex expressions I need to derive from it.
1/sqrt(b) = sqrt(1)/sqrt(b) = sqrt(1/b) is not true when b is negative. Consider b=-1: 1/sqrt(-1) = 1/i = i/(i*i) = i/-1 = -i whereas sqrt(1/-1) = sqrt(-1) = i. If you want this replacement to happen I don't think you have any choice but to do the replacement f.subs(1/sqrt(b), sqrt(1/b)):
>>> f.subs(1/sqrt(b),sqrt(1/b)).subs([(a,1),(b,-1e-9)]).evalf(2)
0.84