integração sympy com função cosseno sob uma raiz quadrada

Estou tentando resolver a integração

integrate( sqrt(1 + cos(2 * x)), (x, 0, pi) )

Claramente, com caneta e papel isso não é difícil, e o resultado é:

Mas ao fazer isso através do Sympy, algo não parece correto.

Tentei os códigos sympy abaixo.

from sympy import *

x = symbols("x", real=True)
integrate(sqrt(1 + cos(2 * x)), (x, 0, pi)).doit()

Então me dá um ValueErrordizendo algo no domínio complexo não definido. Mas eu já defini o símbolo x como uma variável no domínio real.

Aqui está o erro completo:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[7], line 4
      1 from sympy import *
      3 x = symbols("x", real=True)
----> 4 integrate(sqrt(1 + cos(2 * x)), (x, 0, pi)).doit()

File C:\Dev_Tools\Anaconda3\Lib\site-packages\sympy\integrals\integrals.py:1567, in integrate(meijerg, conds, risch, heurisch, manual, *args, **kwargs)
   1564 integral = Integral(*args, **kwargs)
   1566 if isinstance(integral, Integral):
-> 1567     return integral.doit(**doit_flags)
   1568 else:
   1569     new_args = [a.doit(**doit_flags) if isinstance(a, Integral) else a
   1570         for a in integral.args]

File C:\Dev_Tools\Anaconda3\Lib\site-packages\sympy\integrals\integrals.py:499, in Integral.doit(self, **hints)
    497 if reps:
    498     undo = {v: k for k, v in reps.items()}
--> 499     did = self.xreplace(reps).doit(**hints)
    500     if isinstance(did, tuple):  # when separate=True
    501         did = tuple([i.xreplace(undo) for i in did])

File C:\Dev_Tools\Anaconda3\Lib\site-packages\sympy\integrals\integrals.py:710, in Integral.doit(self, **hints)
    707 uneval = Add(*[eval_factored(f, x, a, b)
    708                for f in integrals])
    709 try:
--> 710     evalued = Add(*others)._eval_interval(x, a, b)
    711     evalued_pw = piecewise_fold(Add(*piecewises))._eval_interval(x, a, b)
    712     function = uneval + evalued + evalued_pw

File C:\Dev_Tools\Anaconda3\Lib\site-packages\sympy\core\expr.py:956, in Expr._eval_interval(self, x, a, b)
    953     domain = Interval(b, a)
    954 # check the singularities of self within the interval
    955 # if singularities is a ConditionSet (not iterable), catch the exception and pass
--> 956 singularities = solveset(self.cancel().as_numer_denom()[1], x,
    957     domain=domain)
    958 for logterm in self.atoms(log):
    959     singularities = singularities | solveset(logterm.args[0], x,
    960         domain=domain)

File C:\Dev_Tools\Anaconda3\Lib\site-packages\sympy\solvers\solveset.py:2252, in solveset(f, symbol, domain)
   2250 if symbol not in _rc:
   2251     x = _rc[0] if domain.is_subset(S.Reals) else _rc[1]
-> 2252     rv = solveset(f.xreplace({symbol: x}), x, domain)
   2253     # try to use the original symbol if possible
   2254     try:

File C:\Dev_Tools\Anaconda3\Lib\site-packages\sympy\solvers\solveset.py:2276, in solveset(f, symbol, domain)
   2273     f = f.xreplace({d: e})
   2274 f = piecewise_fold(f)
-> 2276 return _solveset(f, symbol, domain, _check=True)

File C:\Dev_Tools\Anaconda3\Lib\site-packages\sympy\solvers\solveset.py:1060, in _solveset(f, symbol, domain, _check)
   1057     result = Union(*[solver(m, symbol) for m in f.args])
   1058 elif _is_function_class_equation(TrigonometricFunction, f, symbol) or \
   1059         _is_function_class_equation(HyperbolicFunction, f, symbol):
-> 1060     result = _solve_trig(f, symbol, domain)
   1061 elif isinstance(f, arg):
   1062     a = f.args[0]

File C:\Dev_Tools\Anaconda3\Lib\site-packages\sympy\solvers\solveset.py:612, in _solve_trig(f, symbol, domain)
    610 sol = None
    611 try:
--> 612     sol = _solve_trig1(f, symbol, domain)
    613 except _SolveTrig1Error:
    614     try:

File C:\Dev_Tools\Anaconda3\Lib\site-packages\sympy\solvers\solveset.py:688, in _solve_trig1(f, symbol, domain)
    685 if g.has(x) or h.has(x):
    686     raise _SolveTrig1Error("change of variable not possible")
--> 688 solns = solveset_complex(g, y) - solveset_complex(h, y)
    689 if isinstance(solns, ConditionSet):
    690     raise _SolveTrig1Error("polynomial has ConditionSet solution")

File C:\Dev_Tools\Anaconda3\Lib\site-packages\sympy\solvers\solveset.py:2284, in solveset_complex(f, symbol)
   2283 def solveset_complex(f, symbol):
-> 2284     return solveset(f, symbol, S.Complexes)

File C:\Dev_Tools\Anaconda3\Lib\site-packages\sympy\solvers\solveset.py:2252, in solveset(f, symbol, domain)
   2250 if symbol not in _rc:
   2251     x = _rc[0] if domain.is_subset(S.Reals) else _rc[1]
-> 2252     rv = solveset(f.xreplace({symbol: x}), x, domain)
   2253     # try to use the original symbol if possible
   2254     try:

File C:\Dev_Tools\Anaconda3\Lib\site-packages\sympy\solvers\solveset.py:2276, in solveset(f, symbol, domain)
   2273     f = f.xreplace({d: e})
   2274 f = piecewise_fold(f)
-> 2276 return _solveset(f, symbol, domain, _check=True)

File C:\Dev_Tools\Anaconda3\Lib\site-packages\sympy\solvers\solveset.py:1110, in _solveset(f, symbol, domain, _check)
   1106     result += _solve_radical(equation, u,
   1107                              symbol,
   1108                              solver)
   1109 elif equation.has(Abs):
-> 1110     result += _solve_abs(f, symbol, domain)
   1111 else:
   1112     result_rational = _solve_as_rational(equation, symbol, domain)

File C:\Dev_Tools\Anaconda3\Lib\site-packages\sympy\solvers\solveset.py:918, in _solve_abs(f, symbol, domain)
    916 """ Helper function to solve equation involving absolute value function """
    917 if not domain.is_subset(S.Reals):
--> 918     raise ValueError(filldedent('''
    919         Absolute values cannot be inverted in the
    920         complex domain.'''))
    921 p, q, r = Wild('p'), Wild('q'), Wild('r')
    922 pattern_match = f.match(p*Abs(q) + r) or {}

ValueError: 
Absolute values cannot be inverted in the complex domain.

Como faço para integrar isso corretamente usando o Sympy?