# Mismatch between numerical and analytic evaluation of an integral

I evaluated the following integral

$$\int_0^1 \sqrt{r} \left | \cos \left(\left(k+\frac{1}{2}\right) \pi r\right)\right | dr$$

NIntegrate[Sqrt[r] Abs[Cos[(k + 1/2) Pi r]], {r, 0, 1}]


getting as a result 0.413232 for k=1

Then I evaluated the same integral as

Assuming[k ∈ Integers && k ⩾ 0,
Integrate[ 1/Sqrt[2] Sqrt[r (1 + Cos[2 (k + 1/2) Pi r])], {r, 0, 1}]]


and then I evaluated the result

-((4 (-1 + FresnelS[1]))/((1 + 2 k)^(3/2) π))


for k=1 and I get the following result 0.137646.

I am checking this result because I am mainly interested in

$$\int_R^1 {\frac{1}{\sqrt2}}\sqrt{r \cos \left(2\left(k+\frac{1}{2}\right) \pi r\right)}dr$$

Integrate[ 1/Sqrt[2] Sqrt[r (1 + Cos[2 (k + 1/2) Pi r])], {r, R, 1}]


as a function of R.

Where is my error?

• does this help int[r_] = Integrate[1/Sqrt[2] Sqrt[r (1 + Cos[2 (k + 1/2) Pi r])], r]; int[R] - int[1] – Algohi Aug 8 '14 at 16:29
• – Michael E2 Aug 9 '14 at 0:08
• It's a good question. It leads to tricky case studies in MMA. I'll come back soon. – Dr. Wolfgang Hintze Aug 10 '14 at 22:08
• Wow, you've been a member for over 2 years and you have not voted once. Seems like you don't get much help from this forum. – RunnyKine Aug 14 '14 at 1:47

After a lengthy study (I'm using version 8) I conclude that there is a bug in Mathematica in the Integrate function when applied to a Sqrt integrand.

Ok. let's go (some patience is required because of the long text)

Let us define the functions corresponding to your integrals.

Remark: because of the relation $1 + cos(2x) = 2 cos^2(x)$ the two forms of the one-parametric integral are identical. But this fact is not used here.

1) Integral with Abs (designated by f)

Here we define four different functions, the numerical integration, one without any assumptions, and two others with the assumptions k∈Integers and k >=0 respectively:

1.1) Numerical integration

f[k_] := NIntegrate[Sqrt[r] Abs[Cos[(k + 1/2) π r]], {r, 0, 1}]


Some values are

Table[{k, f[k]}, {k, 0, 5, 1/2}]

 {{0, 0.357615}, {1/2, 0.413572}, {1, 0.413232}, {3/2, 0.420072}, {2, 0.419477},
{5/2, 0.421924}, {3, 0.421518}, {7/2, 0.422747}, {4, 0.422465}, {9/2, 0.423197},
{5, 0.422991}}


1.2) Integral without Assumptions

This remains unevaluated by Mathematica.

f0[k_] = Integrate[Sqrt[r] Abs[Cos[(k + 1/2) π r]], {r, 0, 1}]

(* suppressing the Output here *)


1.3) Integral with domain Assumptions

These two integrals are evaluated explicitly, and they give the same reult:

f0ai[k_] =
Integrate[Sqrt[r] Abs[Cos[(k + 1/2) π r]], {r, 0, 1},
Assumptions -> k ∈ Integers]

(-2 Sqrt[1 + 2 k] Cot[k π] +
2 Csc[k π] FresnelS[Sqrt[1 + 2 k]])/((1 + 2 k)^(3/2) π Sqrt[
Csc[k π]^2])

f0ar[k_] =
Integrate[Sqrt[r] Abs[Cos[(k + 1/2) π r]], {r, 0, 1},
Assumptions -> k ∈ Reals]

(-2 Sqrt[1 + 2 k] Cot[k π] +
2 Csc[k π] FresnelS[Sqrt[1 + 2 k]])/((1 + 2 k)^(3/2) π Sqrt[
Csc[k π]^2])


In the domain Complexes the integral remains unevaluated:

f0ac[k_] =
Integrate[Sqrt[r] Abs[Cos[(k + 1/2) π r]], {r, 0, 1},
Assumptions -> k ∈ Complexes]

Integrate[Sqrt[r] Abs[Cos[(1/2 + k) π r]], {r, 0, 1},
Assumptions -> k ∈ Complexes]


1.4) Integral with range Assumptions

These are evaluated explicitly but they are different from those with a domain Assumption

f0age0[k_] =
Integrate[Sqrt[r] Abs[Cos[(k + 1/2) π r]], {r, 0, 1},
Assumptions -> k >= 0]

(4 - 4 FresnelS[1] +
2 Sqrt[Csc[k π]^2] FresnelS[Sqrt[1 + 2 k]] Sin[k π] -
Sqrt[(1 + 2 k) Csc[k π]^2] Sin[2 k π])/((1 + 2 k)^(3/2) π)

f0agt0[k_] =
Integrate[Sqrt[r] Abs[Cos[(k + 1/2) π r]], {r, 0, 1},
Assumptions -> k > 0]

(4 - 4 FresnelS[1] +
2 Sqrt[Csc[k π]^2] FresnelS[Sqrt[1 + 2 k]] Sin[k π] -
Sqrt[(1 + 2 k) Csc[k π]^2] Sin[2 k π])/((1 + 2 k)^(3/2) π)


The integral with k<0 one is different from those with k>=0 or k>0:

f0alt0[k_] =
Integrate[Sqrt[r] Abs[Cos[(k + 1/2) π r]], {r, 0, 1},
Assumptions -> k < 0]

(-2 Sqrt[1 + 2 k] Cot[k π] +
2 Csc[k π] FresnelS[Sqrt[1 + 2 k]])/((1 + 2 k)^(3/2) π Sqrt[
Csc[k π]^2])


Plotting all functions defined so far

Plot[{g[k], f[k], f0ai[k], f0ar[k], f0age0[k], f0alt0[k]}, {k, -6, 5},
PlotRange -> {-1, 1}]


(* plot fs *)

Summarizing up to this point we have

(i) Assumptions -> k ∈ Integers or k ∈ Integers leads to explicit results which are, however, obviously wrong because they are negative in some regions of k while the integrand is definitely positive.

(ii) The Assumption k>=0 leads also to an explicit result which is, however, different from that with the domain Assumptions but still it is wrong because it can become negative. But it seems that the part for 0<=k<1 is ok.

(iii) for k < 0 the functions with are explicit expressions behave "slightly better" but still they are wrong.

In my opinion we have found several bugs in the function Integrate in Mathematica (version 8).

2) Integral with Sqrt[1+Cos] (designated by g)

g[k_] := NIntegrate[1/Sqrt[2] Sqrt[r (1 + Cos[2 (k + 1/2) Pi r])], {r, 0, 1}]

g0[k_] = Integrate[1/Sqrt[2] Sqrt[r (1 + Cos[2 (k + 1/2) Pi r])], {r, 0, 1}]

-((2 (Sqrt[1 + 2 k] Cot[k π] - Csc[k π] FresnelS[Sqrt[1 + 2 k]]) Sqrt[
Sin[k π]^2])/((1 + 2 k)^(3/2) π))

Plot[{g[k], g0[k]}, {k, -5, 5}, PlotRange -> {-1, 1}]


(* suppressing error Messages *)

(* Image plot gs *)

Discussion

(i) The NIntegrate function g[k] is identical to f[k] (ii) the explicit integral g0[k] suffers from the same bugs as discussed for f

3) Integral with Sqrt[1+Cos] and variable lower integration limit R (designated by h)

h[k_, R_] :=
NIntegrate[1/Sqrt[2] Sqrt[r (1 + Cos[2 (k + 1/2) Pi r])], {r, R, 1}]

Plot3D[h[k, R], {R, 0, 1}, {k, -5, 5}]


(* suppressing error Messages ) ( plot h *)

h0[k_, R_] =
Integrate[1/Sqrt[2] Sqrt[r (1 + Cos[2 (k + 1/2) Pi r])], {r, R, 1}]


$Aborted Not evaluated in reasonable time. 4) Simplification of bug discussion We can study a more simplified version of the integral which exhibits the same features: u[k_] := NIntegrate[Sqrt[1 + Cos[2 π k r]], {r, 0, 1}] u0[k_] = Integrate[Sqrt[1 + Cos[2 π k r]], {r, 0, 1}] (Sqrt[1 + Cos[2 k π]] Tan[k π])/(k π) Plot[{u[k], u0[k]}, {k, -5, 5}, PlotRange -> {-1, 2}]  (* plot u *) The problem comes from the Sqrt in the integrand. If we remove it the problem vanishes as can be seen thus: v[k_] := NIntegrate[(1 + Cos[2 π k r]), {r, 0, 1}] v0[k_] = Integrate[(1 + Cos[2 π k r]), {r, 0, 1}] 1 + Sin[2 k π]/(2 k π) Plot[{v[k], v0[k] - 0.1}, {k, -5, 5}, PlotRange -> {0, 2}]  (* plot vs *) The function Sqrt[f[z]] has branch points in the complex z plane at the position of the zeroes of the function f[z]. In our example these are located at Cos[2π k r] == -1 leading to 2π k r = π(n+1/2) or r = (1/2k)(n+1/2)... But I couldn't identify any cause for the buggy behaviour. The discussion resembles one with Daniel Lichtblau some time ago about Integrate in the complex plane. Maybe Daniel can help. Regards, Wolfgang • I have found the cause of the mismatch, and I have found a remedy. As a result, the behaviour of MMA need not be considered a bug, but the function Integrate should be improved, and I'll soon propose how. – Dr. Wolfgang Hintze Aug 11 '14 at 17:58 • Dr. Hintze, it seems that behavior is a bit different in version 10.1. I hope that you will update this Q&A with new information when you have time. Thanks again for your contributions. – Mr.Wizard May 26 '15 at 12:43 In this second answer I give the cause for the mismatch in the integrals, show how to remove it, and make a suggestion to improve the function Integrate[]. Simplified restatement of the problem In order to focus on the core of the problem we consider the simpler integral$\int_0^1 \sqrt{\cos (2 π k r)+1} \, dr$. It has the square root and the cosine with the apropiate argument. Let us define some functions related to the integral in question: 1) Numerical Integration un[k_] := NIntegrate[Sqrt[1 + Cos[2 π k r]], {r, 0, 1}]  k must be asigned a numeric value. Otherwise an error message pops up: un[k]  During evaluation of In[230]:= NIntegrate::inumr: The integrand Sqrt[1+Cos[2 k π r]] has evaluated to non-numerical values for all sampling points in the region with boundaries {{0,1}}. >> un[1] 0.900316  2) Delayed assignment ud[k_] := Integrate[Sqrt[1 + Cos[2 π k r]], {r, 0, 1}]  Case 1: k is left a symbol The integral is calculated analytically. The result is the same as for the immediate assignment ux[k]. ud[k] (Sqrt[1 + Cos[2 k π]]] Tan[k π])/(k π)  Case 2: k has a numeric value The integral is calculated analytically with the given value of k ud[1] (2 Sqrt[2])/π  3) Immediate assignment ux[k_] = Integrate[Sqrt[1 + Cos[2 π k r]], {r, 0, 1}] (Sqrt[1 + Cos[2 k π]] Tan[k π])/(k π)  Case 1: k is left symbolic (see above) Case 2: k has a numeric value (see below) ux[1/3] (3 Sqrt[3/2])/π  Now we calculate the values at k=1 for the three functions: un[1] 0.900316 ud[1] % // N (2 Sqrt[2])/π 0.900316 ux[1] 0  Now, for comparison, we plot the three functions as a function of k: ad 1) Numerical integration (this is the "true" result of the integral) Plot[{un[k]}, {k, 0, 5}, PlotRange -> {-1, 2}, PlotLabel -> "un[k]", AxesLabel -> {"k", "un[k]"}] (* plot_un * )  ad2) Delayed assignment Leads to correct results but requires exceedingly long time. Therefore we have cheated a bit using ListLinePlot insted of Plot ListLinePlot[Table[{k, ud[k]}, {k, 0, 5, 1/10}], PlotRange -> {-1, 2}, PlotLabel -> "ud[k]", AxesLabel -> {"k", "ud[k]"}] (* plot_ud *)  ad 3) Immeadiate assignment This is a plot of the explicit analytic calculation of the integral with symbolic k Plot[{ux[k]}, {k, 0, 5}, PlotRange -> {-1, 2}, PlotLabel -> "ux[k]", AxesLabel -> {"k", "ux[k]"}] (* plot_ux *)  We can clearly see the "mismatch" between the true solutions (1 and 2) and the analytic expression with symbolic k (3). This is the problem to be explained and resolved. Study of the mismatch In order to understand the mismatch we need to look at the way a definite integral is calculated in Mathematica. It is done in two steps step 1: calculate the indefinite integral step 2: insert the start and end value of the integration variable into the result of step 1 and subtract the first quantity from the second Some comments are apropriate here ad step 1 The indefinite integral is the so called antiderivative (Stammfunktion) of the integrand. This is an expression which on differentiating gives back the integrand. It is clear that the antiderivative is defined only up to an additive constant expression (with respect to the integration variable). Hence it is better to speak of "an antiderivative". We shall call the antiderivative provided by Integrate the "naive antiderivative". ad step 2 This is a simplified description. We must be very careful here. Because here lies the cause of a possible mismatch. The common procedure described abowe in step 2 is only correct under the following condition: the antidrivative must be continuous along the integration path. This condition is normally not mentioned in elementary math textbooks. But here we have found an example where the "naive" antiderivative obtained by Integrate is not continous, and wrong results pop up. The "naive" antiderivative in our example is given by ad0[r_, k_] = Integrate[Sqrt[1 + Cos[2 π k r]], r] (Sqrt[1 + Cos[2 k π r]] Tan[k π r])/(k π)  We need to investigate the continuity with respect to the integration variable for given parameter k. For k = 1 we have ad0[r, 1] (Sqrt[1 + Cos[2 π r]] Tan[π r])/π Plot[ad0[r, 1], {r, 0, 1}, PlotRange -> {-1, 2}, PlotLabel -> "ad0: 'naive' antiderivative\nas a function of the integration \ variable r for k=1", AxesLabel -> {"r", "ad0[r,k=1]"}] (* plot_ad0 *)  Obviously the antiderivative has a discontinuity at r=1/2. This is due to the joint effect of two singularities at this point: a pole from Tan and a branch point from the Sqrt. The size d of the jump is the difference between the limits from below and above 1/2: Limit[ad0[r, 1], r -> 1/2, Direction -> +1] Sqrt[2]/π Limit[ad0[r, 1], r -> 1/2, Direction -> -1] -(Sqrt[2]/π) d = %% - % (2 Sqrt[2])/π  By adding this jump we can obtain an antiderivative which is continuous at r=1/2. This is done by defining (more generally, not fixing k to 1) ad1[r_, k_] := ad0[r, k] /; r k <= 1/2 ad1[r_, k_] := ad0[r, k] + 2 Sqrt[2]/(k π) /; 1/2 < r k <= 3/2  Plotting it shows the continuity Plot[ad1[r, 1], {r, 0, 1}, PlotRange -> {-1, 2}, PlotLabel -> "ad1: antiderivative with one jump removed\nas a function of the \ integration variable r for k=1", AxesLabel -> {"r","ad1[r,k=1]"}] (* plot_ad1r *)  We can now apply the simple rule of inserting end and start values into the antiderivative and subtracting to get the values of the integral from r=0 to r=1 for k<3/2 ad1[1, 1] - ad1[0, 1] (2 Sqrt[2])/π  This value is correct as can be show by comparison with ud1: ud[1] (2 Sqrt[2])/π  And also over the range of k covered so far (0<k<3/2): Plot[ad1[1, k], {k, 0, 5}, PlotRange -> {-1, 2}, PlotLabel -> "ad1: antiderivative with one jump removed\nas a function of k with the \ integration variable r=1", AxesLabel -> {"k","ad1[r=1,k]"}] (* plot_ad1k *)  This procedure can be generalized to generate a smooth antiderivative up to a given value of r*k. Indeed using ad2[r_, k_] := Piecewise[Table[{ad0[r, k] + m Sqrt[2]/(k π), (m - 1)/2 < r k <= (m + 1)/2}, {m, 0, 8, 2}]]  we have defined ad2[r,k] $\begin{array}{ll} \{ & \begin{array}{ll} \frac{\sqrt{1+\text{Cos}[2 k \pi r]} \text{Tan}[k \pi r]}{k \pi } & -\frac{1}{2}<k r\leq \frac{1}{2} \\ \frac{2 \sqrt{2}}{k \pi }+\frac{\sqrt{1+\text{Cos}[2 k \pi r]} \text{Tan}[k \pi r]}{k \pi } & \frac{1}{2}<k r\leq \frac{3}{2} \\ \frac{4 \sqrt{2}}{k \pi }+\frac{\sqrt{1+\text{Cos}[2 k \pi r]} \text{Tan}[k \pi r]}{k \pi } & \frac{3}{2}<k r\leq \frac{5}{2} \\ \frac{6 \sqrt{2}}{k \pi }+\frac{\sqrt{1+\text{Cos}[2 k \pi r]} \text{Tan}[k \pi r]}{k \pi } & \frac{5}{2}<k r\leq \frac{7}{2} \\ \frac{8 \sqrt{2}}{k \pi }+\frac{\sqrt{1+\text{Cos}[2 k \pi r]} \text{Tan}[k \pi r]}{k \pi } & \frac{7}{2}<k r\leq \frac{9}{2} \\ 0 & \text{True} \\ \end{array} \\ \end{array} $Showing the antiderivatives in 3D plots Plot3D[ad0[r, k], {r, 0, 1}, {k, 0, 5}, PlotRange -> {-1, 1.5}, PlotLabel -> "'naive' antiderivative: discontinious at r k = n+1/2, n>=1", AxesLabel -> {"r", "k"}] (* plot3d_ad0 *)  Plot3D[ad1[r, k], {r, 0, 1}, {k, 0, 5}, PlotRange -> {-1, 1.5}, PlotLabel -> "antiderivative: two discontinuities removed up to r k <= 3/2", AxesLabel -> {"r", "k"}] (* plot3d_ad1 *)  Plot3D[ad2[r, k], {r, 0, 1}, {k, 0, 5}, PlotRange -> {-1, 1.5}, PlotLabel -> "antiderivative: four discontinuities removed up to r k <= 9/2\nThe \ white stripes are an artefact from Piecewise", AxesLabel -> {"r", "k"}] (* plot3d_ad2 *)  Summary We have shown the cause of discrepancies of values of a specific integral$\int_0^1 \sqrt{\cos (2 \pi k r)+1} \, dr$depending on the usage of the function Integrate[]. This discrepancy arises when the "naive" antiderivative, i.e. the analytic result of Integrate[], is not continous on the path of integration. We have also shown (in this example) how to generate an antiderivative which is smooth for a certain region of the parameters. These results generalise and can be applied to the original integrals of Steffano's problem. I would not consider this behaviour a bug in Mathematica, but it would be very helpful if the function Integrate could be improved to give a warning when discontinous antiderivatives appear. As to Steffanos's question: your error was in the interpretation of some of the the commands you were using. In a sense you stepped into a trap provided by Mathematica. But I am grateful for your question. It was fun to study the problem, and I hope it helps you in your work. Best regards, Wolfgang Addendum 12.08.14 21:20: exact solution of the integral in terms of elliptic integrals Studying possible regularizations I discovered that our integral$\int _0^1\sqrt{\cos (2 \pi k r)+1}dr$can be expressed exactly by the elliptic integral of the second kind fi[k_] := Sqrt[2] EllipticE[π k, 1]/(π k)  This function is valid for all values of k. The former exakt solution ux[k] = (Sqrt[1 + Cos[2 k π]] Tan[k π])/(k π) coincides with it only in the interval -1/2<k<1/2 Let us compare both in a plot Plot[{fi[k], ux[k]}, {k, 0, 5}, PlotRange -> {-1, 2}, PlotLabel -> "fi[k]: exact analytic solution of the integral in terms of a \n\ complete elliptic integral of the second kind, valid for all k.\nAnd \ ux[k], the former exact solution, which is valid only for |k|<1/2", AxesLabel -> {"k", "fi[k], ux[k]"}] (* plot_fi _ux *)  Furthermore the function adf[r_, k_] := Sqrt[2 ] EllipticE[π k r, 1]/(π k)  is the antiderivative which is continuous in r for all values of k. In fact, differentiating with respect to r gives D[Sqrt[2] EllipticE[k r π, 1]/(π k), r] // Simplify Sqrt[1 + Cos[2 k π r]]  Funny, the elliptic integral function incorporates the rather complicated piecewiese removal of jumps. To explore the relation between ux and fi seems to be an interesting pasttime. It is not analytic continuation but removal of discontinuities which relates them ... At least it is fair to say that the elliptic integral is a much better analytic solution than ux. And MMA should tell us the limits of validity of ux[k] as the exact solution of the integral. PS 13.08.14 The integral can be found (of course) in Gradsteyn/Ryshik (5th Edition) 2.614.2: E(x,1) Regards, Wolfgang • A brief look seems to indicate that my analysis explains some of the related Integrate[]-questions here in StackExchange as well. I need to go more deeply into the subject, of course. – Dr. Wolfgang Hintze Aug 12 '14 at 16:08 • +1 extremely instructive and fun reading. How did you discover the EllipticE-alternative? Hopefully not with MMA :) – eldo Aug 12 '14 at 20:02 • @ eldo: thanks, it was a pleasure. I found the EllipticE alternative by regularizing the integral thus u[k,a] = Integrate[Sqrt[1+a Cos[2 Pi k r],{r,0,1}]. Here MMA Returns the E-function, and there is no Problem letting a->1. Hence it was infact MMA (with some help from my side) ;) – Dr. Wolfgang Hintze Aug 13 '14 at 6:50 • I am pleased you chose this site as your new "home" for such insights. Thank you! – Mr.Wizard Aug 15 '14 at 0:07 • Thank you I solved the problem few days ago using the Bessel form of the integrand function. Best Regards – Raffaele Carlone Aug 15 '14 at 17:43 The exact analytic soultion 1. Introduction The problem was still intriguing me with the result of a further study which I present in the following, for clarity as another solution. I have chosen to write the formulas in the more theoretical text in traditonal form. Abstract We calculate here the explcit analytic solution for the integral$f(k,R)=\int _R^1r^a |\cos (\pi r k)|dr$with$a=0$or$a=\frac{1}{2}$, respectively, and$0\leq R<1$. This covers the two integrals in the orginal formulation of the problem. The results for$R=0$are$f(k, a=0)=\frac{2 \left\lfloor k+\frac{1}{2}\right\rfloor +(-1)^{\left\lfloor k+\frac{1}{2}\right\rfloor } \sin (\pi k)}{\pi k}f(k, a=1/2)=\frac{\sum _{m=1}^{\left\lfloor k+\frac{1}{2}\right\rfloor } \left(\sqrt{2} (-1)^m S\left(\sqrt{2 m-1}\right)+\sqrt{4 m-2}\right)+(-1)^{\left\lfloor k+\frac{1}{2}\right\rfloor } \left(\sqrt{k} \sin (\pi k)-\frac{S\left(\sqrt{2} \sqrt{k}\right)}{\sqrt{2}}\right)}{\pi k^{3/2}}$Here$S()$is the function FresnelS. The results for$0\leq R<1$can be expressed as$f(k,R)=f(k)-\left(\frac{R}{k}\right)^{a+1} f(k R)$2. Main part Changing the integration variable to$t\to \pi k r$leads to$f(k)=\frac{g(k)}{(\pi k)^{a+1}}$where$g(k)=\int _{0}^{\pi k}t^a |\cos (t)|dt$This can be made more explicit by considering the intervals of k where cos > 0 and cos < 0 for increasing k. We shall call thees functions valid in certain interval in k "partialwaves". Here are the first few$g(0;k)=\int_0^{\pi k} t^a \cos (t) \, dt,0\leq k<\frac{1}{2}g(1;k)=\int_0^{\frac{\pi }{2}} t^a \cos (t) \, dt-\int_{\frac{\pi }{2}}^{\pi k} t^a \cos (t) \, dt,\frac{1}{2}\leq k<\frac{3}{2}g(2;k)=\int _0^{\frac{\pi }{2}}t^a \cos (t)dt-\int _{\frac{\pi }{2}}^{\frac{3 \pi }{2}}t^a \cos (t)dt + \int _{\frac{3 \pi }{2}}^{\pi k}t^a \cos (t)dt,\frac{3}{2}\leq k<\frac{5}{2}$and so on. This can be written recursively for n = 1, 2, 3, ... as$g(n;k)=g\left(n-1;n-\frac{1}{2}\right)+(-1)^n \int _{\pi \left(n-\frac{1}{2}\right)}^{\pi k}t^a \cos (t)dt,n-\frac{1}{2}\leq k<n+\frac{1}{2}$In order to calculate the functions g(n;k) (designated ggX[n,k], where X->0 for a = 0 and X->1 for a = 1/2) we need first the expressions for n=0: Case a = 0 gg0[0, k_] = Integrate[Cos[t], {t, 0, π k}] Sin[k π ]  Case a = 1/2 gg1[0, k_] = Integrate[Sqrt[t] Cos[t], {t, 0, π k}] -Sqrt[(π /2)] FresnelS[Sqrt[2] Sqrt[k]] + Sqrt[k] Sqrt[π] Sin[k π ]  Now for n = 1, 2, ... the part of g denpending on k ist the last integral which, including the factor (-1)^n, will be designated gfX[n,k] For a = 0 we find Simplify[(-1)^n Integrate[t^a Cos[t], {t, π (n - 1/2), π k}, Assumptions -> { a == 0}], {n \[Element] Integers }] (-1)^n ((-1)^n + Sin[k π])  (helping Mathematica a bit in the trivial simplification of the minus ones) gf0[n_, k_] = 1 + (-1)^n Sin[k π] 1 + (-1)^n Sin[k π]  and for n > 0 and a = 1/2 gf1[n_, k_] = Simplify[(-1)^n Integrate[ Sqrt[t] Cos[t], {t, π (n - 1/2), π k}], {n \[Element] Integers , n > 0, k > n - 1/2}] 1/2 (-1)^n Sqrt[π] ((-1)^n Sqrt[-2 + 4 n] - Sqrt[2] FresnelS[Sqrt[2] Sqrt[k]] + Sqrt[2] FresnelS[Sqrt[-1 + 2 n]] + 2 Sqrt[k] Sin[k π])  The recursion relation is for a = 0 Clear[gg0] req = {gg0[0, k] == Sin[k π], gg0[n, k] == gg0[n - 1, π (n - 1/2)] + gf0[n, k]} {gg0[0, k] == Sin[k π], gg0[n, k] == 1 + gg0[-1 + n, (-(1/2) + n) π] + (-1)^n Sin[k π]}  The direct attack RSolve[req, gg0[n, k], n] RSolve[{gg0[0, k] == Sin[k π], gg0[n, k] == 1 + gg0[-1 + n, (-(1/2) + n) π] + (-1)^n Sin[k π]}, gg0[n, k], n]  does not work. Therefore we do it with a the inductive method "examples and guess" as follows: gg0[0, k_] = Sin[k π]; Table[x = gg0[n - 1, (n - 1/2)] + gf0[n, k]; gg0[n_, k_] = x, {n, 1, 5}] {2 - Sin[k π], 4 + Sin[k π], 6 - Sin[k π], 8 + Sin[k π], 10 - Sin[k π]}  The general formula is easily deduced to be Clear[gg0] gg0[n_, k_] = (2 n + (-1)^n Sin[π k]);  which holds even for n=0. The "partialwaves" ff0[n,k] are therefore given by ff0[n_, k_] = gg0[n, k]/(π k) (2 n + (-1)^n Sin[k π])/(k π)  It is interesting to plot some of them together in one diagram Plot[Evaluate[Table[ff0[n, k], {n, 0, 3}]], {k, 0, 7/2}, PlotRange -> {-0.5, 1.5}, PlotLabel -> "'Partialwaves'(0..3) of \!$$\*TagBox[\(f0(k) = \ \*SubsuperscriptBox[\(\[Integral]$$, $$0$$, \ $$1$$]\*TemplateBox[{RowBox[{\"cos\", \"(\", \nRowBox[{\"\[Pi]\", \" \ \", \"k\", \" \", \"r\"}], \")\"}]},\n\"Abs\"] \[DifferentialD]r\), HoldForm]\)", AxesLabel -> {"k", "f0(k)"} ]  (* plot_ff0_pw *) We see that the functions neatly play the "relay"-game in defining the solution. The blue curve, valid for 0<=k<1/2, relays to the red one at k=1/2. The red has the defining lead from k=1/2 to k=3/2 where it finally relays to the yellow curve, and so forth. Now we take the final step relating n to k via the obvious formula nn[k_] := Floor[k + 1/2]  The complete analytic solution of the initial integral is therefore for a = 0 f0[k_] = ff0[Floor[k + 1/2], k];  as mentioned in the beginning. The asymptotic behaviour of f0 for k->[Infinity] is obviously f0as = 2/π % // N 2/π 0.63662  In the case a = 1/2 the inductive method gives (To get simpler expressions we divide by Sqrt[π]) gg1[0, k]/Sqrt[π] // Simplify -(FresnelS[Sqrt[2] Sqrt[k]]/Sqrt[2]) + Sqrt[k] Sin[k π] 1/Sqrt[π] Table[x = gg1[n - 1, (n - 1/2)] + gf1[n, k]; gg1[n_, k_] = x, {n, 1, 4}] // Expand // TableForm $\begin{array}{l} \sqrt{2}-\sqrt{2} \text{FresnelS}[1]+\frac{\text{FresnelS}\left[\sqrt{2} \sqrt{k}\right]}{\sqrt{2}}-\sqrt{k} \text{Sin}[k \pi ] \\ \sqrt{2}+\sqrt{6}-\sqrt{2} \text{FresnelS}[1]+\sqrt{2} \text{FresnelS}\left[\sqrt{3}\right]-\frac{\text{FresnelS}\left[\sqrt{2} \sqrt{k}\right]}{\sqrt{2}}+\sqrt{k} \text{Sin}[k \pi ] \\ \sqrt{2}+\sqrt{6}+\sqrt{10}-\sqrt{2} \text{FresnelS}[1]+\sqrt{2} \text{FresnelS}\left[\sqrt{3}\right]-\sqrt{2} \text{FresnelS}\left[\sqrt{5}\right]+\frac{\text{FresnelS}\left[\sqrt{2} \sqrt{k}\right]}{\sqrt{2}}-\sqrt{k} \text{Sin}[k \pi ] \\ \sqrt{2}+\sqrt{6}+\sqrt{10}+\sqrt{14}-\sqrt{2} \text{FresnelS}[1]+\sqrt{2} \text{FresnelS}\left[\sqrt{3}\right]-\sqrt{2} \text{FresnelS}\left[\sqrt{5}\right]+\sqrt{2} \text{FresnelS}\left[\sqrt{7}\right]-\frac{\text{FresnelS}\left[\sqrt{2} \sqrt{k}\right]}{\sqrt{2}}+\sqrt{k} \text{Sin}[k \pi ] \\ \end{array}$We de deduce the formula Clear[gg1] gg1[n_, k_] := Sum[Sqrt[4 m - 2] + (-1)^m Sqrt[2] FresnelS[Sqrt[2 m - 1]], {m, 1, n}] - (-1)^n (FresnelS[Sqrt[2] Sqrt[k]]/Sqrt[2] - Sqrt[k] Sin[k π])  The functions ff1 are ff1[n_, k_] = gg1[n, k]/(π k^(3/2))  Plotting the first 4 partialwaves: Plot[Evaluate[Table[ff1[n, k], {n, 0, 3}]], {k, 0, 7/2}, PlotRange -> {-0.5, 1.5}, PlotLabel -> "'Partialwaves'(0..3) of \!$$\*TagBox[\(f1 \((k)$$ = \ \*SubsuperscriptBox[$$\[Integral]$$, $$0$$, $$1$$]\*SqrtBox[$$r$$] \ \*TemplateBox[{RowBox[{\"cos\", \"(\", \nRowBox[{\"\[Pi]\", \" \", \ \"k\", \" \", \"r\"}], \")\"}]},\n\"Abs\"] \[DifferentialD]r\), HoldForm]\)", AxesLabel -> {"k", "f1(k)"} ]  (* plot_ff1_pw *) A similar "relay"-game is observed. Notice that the dependence of f on k is given by Sin[π k]/(π k), exactly as in the case a=0. The final explicit analytical solution is obtained as before by replacing n via k. Is has the form a finite sum: f1[k_] := 1/(π k^(3/2)) (Sum[ Sqrt[4 m - 2] + (-1)^m Sqrt[2] FresnelS[Sqrt[2 m - 1]], {m, 1, Floor[k + 1/2]}] - (-1)^ Floor[k + 1/2] (FresnelS[Sqrt[2] Sqrt[k]]/Sqrt[2] - Sqrt[k] Sin[k π])) // Simplify  as mentioned in the beginning. In order to caluclate the asymptotic behaviour of f1[k] for k->[Infinity] we observe that the only term left over is 1/(π k^(3/2)) Sum[Sqrt[4 m - 2], {m, 1, Floor[k + 1/2]}]  The in the limit k->[Infinity] this goes to f1as = 4/(3 π ); % // N 0.424413  The (not so difficult) proof is left as an exercise to the reader. Now we plot the functions together with their respective asymptotic value Plot[{2/π], f0[k]}, {k, 0, 5}, PlotRange -> {0, 1.1}, PlotLabel -> "The function \!$$\*TagBox[\(f0 \((k)$$ = \*SubsuperscriptBox[$$\[Integral]\$$, $$0$$, $$1$$]\*TemplateBox[{RowBox[{\"cos\", \"(\", \nRowBox[{\"\[Pi]\", \ \" \", \"k\", \" \", \"r\"}], \")\"}]},\n\"Abs\"] \[DifferentialD]r\), HoldForm]\)\nand the asymptotic value \!$$\*TagBox[FractionBox[\(2$$, $$\[Pi]\$$], HoldForm]\)", AxesLabel -> {"k", "f0(k)"}]  (* plot_f0 *) Plot[{4/(3 π), f1[k]}, {k, 0, 5}, PlotRange -> {0, 1.1}, PlotLabel -> "The function \!$$\*TagBox[\(f1 \((k)$$ = \*SubsuperscriptBox[$$\ \[Integral]$$, $$0$$, $$1$$]\*SqrtBox[$$r$$] \ \*TemplateBox[{RowBox[{\"cos\", \"(\", \nRowBox[{\"\[Pi]\", \" \", \ \"k\", \" \", \"r\"}], \")\"}]},\n\"Abs\"] \[DifferentialD]r\), HoldForm]\)\nand the asymptotic value \!$$\*TagBox[FractionBox[\(4$$, \ $$3\\\ \[Pi]$$], HoldForm]\)", AxesLabel -> {"k", "f1(k)"}]  (* plot_f1 *) 3. Summary 1) The exact analytic solution of the integral$f(k)=\int _0^1r^a |\cos (\pi r k)|dr$was calculated with the result$f(k, a=0)=\frac{2 \left\lfloor k+\frac{1}{2}\right\rfloor +(-1)^{\left\lfloor k+\frac{1}{2}\right\rfloor } \sin (\pi k)}{\pi k}$It is the sum of two highly dicontinuious functions which merge to give a smooth function. For a=1/2 the Situation is similar. 2) The solutions are oscillating about a positive value. The oscillation is goverened by the well-known function$\frac{\sin (x)}{x}$3) As we have shown earlier, the solution is also given by the elliptic integral, so that we have found here an interesting expression (which I coulnd't find in Abramovich/Stegun)$E(\pi k|1) = 2 \left\lfloor k+\frac{1}{2}\right\rfloor +(-1)^{\left\lfloor k+\frac{1}{2}\right\rfloor } \sin (\pi k)$4) The procedure described here generalizes to other values of the Parameter a > -1. 5) Last but not least let me point out that I couldn't have found the results without MMA - in finite time! This might justify the presentation of the more or less mathematical topic in this forum. • I have the result for 0<R<1. We can generalize a little bit the result. I agree with your solution. – Raffaele Carlone Aug 20 '14 at 15:16 • @ Raffaele: please present your solution here so we can discuss it. Thank you. – Dr. Wolfgang Hintze Aug 20 '14 at 16:26 just an extended comment: In Win7, MMA7 your original question is solved correctly k = 1 r1 = Block[{r}, NIntegrate[Sqrt[r] Abs[Cos[(k + 1/2) Pi r]], {r, 0, 1}]] r2 = Assuming[k \[Element] Integers && k \[GreaterSlantEqual] 0, Integrate[1/Sqrt[2] Sqrt[r (1 + Cos[2 (k + 1/2) Pi r])], {r, 0, 1}]] // N rDiff = r1 - r2 (*out*) 1 0.413231743899931 0.4132317438999164 1.454392162258955*10^-14  ## Answer to comment @hintze The assumptions made in Assuming before and those made by Integrate are different k =. Block[{r}, NIntegrate[Sqrt[r] Abs[Cos[(k + 1/2) Pi r]], {r, 0, 1}]] /. k -> 1. Integrate[ 1/Sqrt[2] Sqrt[r (1 + Cos[2 (k + 1/2) Pi r])], {r, 0, 1}] /. k -> 1. (* out *) NIntegrate::inumr: The integrand Sqrt[r] Abs[Cos[(1/2+k) \[Pi] r]] has evaluated to non-numerical values for all sampling points in the region with boundaries {{0,1}}. >> NIntegrate::inumr: The integrand Sqrt[r] Abs[Cos[(1/2+k) \[Pi] r]] has evaluated to non-numerical values for all sampling points in the region with boundaries {{0,1}}. >> 0.413231743899931 0.4132317438992997  Wolfgang's reply to hieron as of 20.08.14 08:40 I continue in the manner hieron started, although such discussion should probably be done using comments, but it is more convenient this way. In $Version

"8.0 for Microsoft Windows (64-bit) (October 7, 2011)"


your second integral is returned unevaluated:

q1 = Integrate[1/Sqrt[2] Sqrt[r (1 + Cos[2 (k + 1/2) Pi r])], {r, 0, 1}] /.
k -> 1

During evaluation of In[1]:= Infinity::indet: Indeterminate expression ComplexInfinity+ComplexInfinity encountered. >>
Indeterminate


We rather need to take the limit, and we do it with the option Direction

q2 = Limit[
Integrate[1/Sqrt[2] Sqrt[r (1 + Cos[2 (k + 1/2) Pi r])], {r, 0, 1}], k -> 1,
Direction -> +1]
% // N

(2 (3 + Sqrt[3] FresnelS[Sqrt[3]]))/(9 \[Pi])

0.275586

q3 = Limit[
Integrate[1/Sqrt[2] Sqrt[r (1 + Cos[2 (k + 1/2) Pi r])], {r, 0, 1}], k -> 1,
Direction -> -1]
% // N

-((2 (3 + Sqrt[3] FresnelS[Sqrt[3]]))/(9 \[Pi]))

-0.275586


Hence we find q3 = - q2.

Now, a similar effect happens at k = 0

q4 = Limit[
Integrate[1/Sqrt[2] Sqrt[r (1 + Cos[2 (k + 1/2) Pi r])], {r, 0, 1}], k -> 0,
Direction -> +1]
% // N

(2 - 2 FresnelS[1])/\[Pi]

0.357615

q5 = Limit[
Integrate[1/Sqrt[2] Sqrt[r (1 + Cos[2 (k + 1/2) Pi r])], {r, 0, 1}], k -> 0,
Direction -> -1]
% // N

(2 (-1 + FresnelS[1]))/\[Pi]

-0.357615


@hieron: please provide the results of these commands in your version. Best regards, Wolfgang

• Ok. But the problem arises if you first calculate the integral leaving k as a symbol and then let k->1. Try it. – Dr. Wolfgang Hintze Aug 19 '14 at 15:16

This is my solution:

n[R_, k_] := k - IntegerPart[(1 - R)/(2/(2 k + 1))];
sum[k_, R_] :=
N[((-1)^n[R, k] (t[k, r[n[R, k], k]] - t[k, R]) +
Sum[(-1)^j (t[k, r[j, k]] - t[k, r[j - 1, k]]), {j, n[R, k] + 1,
k}])];
r[x_, k_] := 1/(2 k + 1) + (2 x )/(2 k + 1);
t[k_, r_] = -((
2 FresnelS[Sqrt[1 + 2 k] Sqrt[r]])/((1 + 2 k)^(3/2) \[Pi])) + (
2 Sqrt[r] Sin[1/2 (1 + 2 k) \[Pi] r])/((1 + 2 k) \[Pi]);

• I have also quite completed the asymptotic estimates, I would ask to Wolfgang Hintze to contact me at my email ozellus@gmail.com. – Raffaele Carlone Aug 23 '14 at 8:51
• @ Raffaele: how did you derive your solution? Also, you might have noticed that I have calculated the asymptotic behaviour. But maybe you wish to discover the results I presented here by yourself. Why not? I have also sent you an email as you requested. – Dr. Wolfgang Hintze Aug 26 '14 at 10:42
• Sorry Wolfang, check the email. – Raffaele Carlone Aug 28 '14 at 11:01