Integrate gives an imaginary part for a 1-d definite integral of a real-valued function

I'm using Integrate to evaluate an integral with known result:

$$\int_0^\pi \frac{\cos(n \theta)-\cos(n \theta_0)}{\cos(\theta)-\cos(\theta_0)} d\theta = \frac{\pi\sin(n \theta_0)}{\sin(\theta_0)},$$

and the code is as follows

int = Assuming[n \[Element] Integers && 0 < t0 < \[Pi] && n > 0,
Integrate[(Cos[n  t] - Cos[n  t0])/(
Cos[t] - Cos[t0]), {t, 0, \[Pi]}]]

(*-(1/2) I \[Pi] Csc[
t0] (2 (Cos[n t0] + (-1)^n Cos[2 n ArcSin[Cos[t0/2]]]) -
I Sqrt[2 \[Pi]] Sqrt[1 - Cos[t0]]
MeijerG[{{1/2}, {1/2 - n, 1/2 + n}}, {{1/2, 1/2}, {0}},
Cos[t0/2]^2])*)

Strangely, it returns a result with an imaginary part. I checked the result with numerical methods, and the result suggests that the real part of the above result is the right result of the integral.

rule = {n -> 3, t0 -> \[Pi]/7};

int /. rule // N

NIntegrate[(Cos[n  t] - Cos[n  t0])/(Cos[t] - Cos[t0]) /. rule //
Evaluate, {t, 0, \[Pi]}]

\[Pi]  Sin[n t0]/Sin[t0] /. rule // N

(*7.05909 - 3.22238 I*)
(*7.05909*)
(*7.05909*)

My problem is: what's wrong with it? How to fix it? My mathematica $Version is "14.0.0 for Mac OS X ARM (64-bit) (December 1, 2023)" Update: With the answer from @Nasser, I found a potential cause. Using integrand = (Cos[n t] - Cos[n t0])/(Cos[t] - Cos[t0]) assume = n ∈ Integers && 0 < t0 < π && n > 0; integrand2 = TrigExpand[integrand] to split the integral into two parts. Both of the parts do not converge, they can be only evaluated in the scene of PrincipalValue -> True. Here we focus on the second part. Using the following code int1 = Integrate[integrand2[[2]], {t, 0, \[Pi]}, Assumptions -> assume, PrincipalValue -> True] int2 = Integrate[integrand2[[2]], {t, 0, \[Pi]}, GenerateConditions -> False, PrincipalValue -> True] Assuming[assume, int1 - int2 // Simplify] (*-I \[Pi] Cos[n t0] Csc[t0]*) (*Cos[n t0] Csc[t0] (Log[-Cot[t0/2]] - Log[Cot[t0/2]])*) (*-2 I \[Pi] Cos[n t0] Csc[t0]*) we can see that the principal values of the integral with different assumptions are different. We expect it to give the result int2 in my assumption, but it gives int1. Maybe it is a bug? • Table[Integrate[(Cos[n t] - Cos[n t0])/(Cos[t] - Cos[t0]), {t, 0, \[Pi]}, Assumptions -> 0 < t0 < \[Pi], GenerateConditions -> False], {n, 1, 10}] performs {\[Pi],2 \[Pi] Cos[t0],\[Pi]+2 \[Pi] Cos[2 t0],\[Pi] Csc[t0] Sin[4 t0],\[Pi] Csc[t0] Sin[5 t0],\[Pi] Csc[t0] Sin[6 t0],\[Pi] Csc[t0] Sin[7 t0],\[Pi] Csc[t0] Sin[8 t0],\[Pi] Csc[t0] Sin[9 t0],\[Pi] Csc[t0] Sin[10 t0]}. Commented Apr 23 at 19:31 2 Answers You seem to have something not quite right with the conditions. Here is a workaround. Split the integrand into 2 parts. Integrate each part on its own and add the result. Use the assumption you had only on the first part. Now you get correct answer. No imaginary part. I noticed this because adding the assumption you had on the second integrand part caused Mathematica to complain it does not converge. integrand = (Cos[n t] - Cos[n t0])/(Cos[t] - Cos[t0]) assume = n ∈ Integers && 0 < t0 < π && n > 0; integrand2 = TrigExpand[integrand]; anti1 = Integrate[integrand2[[1]], {t, 0, π}, Assumptions -> assume]; anti2 = Integrate[integrand2[[2]], {t, 0, π}]; (anti1 + anti2[[1]]) /. {n -> 3, t0 -> π/7} N[%] // Chop Gives (*7.05909*) You might want to report this to Wolfram support. • Both of the parts don't converge, They can be evaluated only in principal value. For example, Mathematica also complains Integrate[integrand2[[2]] /. n -> 2.5 /. t0 -> 1/4, {t, 0, \[Pi]}] does not converge. Commented Apr 23 at 15:27 An alternate approach$Version

(* "14.0.0 for Mac OS X ARM (64-bit) (December 13, 2023)" *)

Clear["Global`*"]

int[n_Integer?Positive] := int[n] =
Assuming[0 < t0 < π,
Integrate[(Cos[n   t] - Cos[n   t0])/(Cos[t] - Cos[t0]),
{t, 0, π}]]

Generate a sequence for several values of n

seq = int /@ Range[6]

(* {π, 2 π Cos[t0], π + 2 π Cos[2 t0],
2 π (Cos[t0] + Cos[3 t0]), π (1 + 2 Cos[2 t0] + 2 Cos[4 t0]),
2 π Cos[t0] (1 + 2 Cos[4 t0])} *)

Use FindSequenceFunction to generalize

int2[n_, t0_] = Assuming[0 < t0 < π,
FindSequenceFunction[seq, n] // FullSimplify]

(* π Csc[t0] Sin[n t0] *)

int2[3, Pi/7] // Simplify

(* 1/2 π Csc[π/14] *)

% // N

(* 7.05909 *)