Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I evaluated the following integral

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

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?

share|improve this question
    
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 at 16:29
    
    
It's a good question. It leads to tricky case studies in MMA. I'll come back soon. –  Dr. Wolfgang Hintze Aug 10 at 22:08
1  
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 at 1:47

5 Answers 5

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}]

enter image description here (* 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 *)

enter image description here (* 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 ) enter image description here ( 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}]

enter image description here

(* 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}]

enter image description here (* 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

share|improve this answer
    
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 at 17:58

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 * )

enter image description here

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 *)

enter image description here

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 *)

enter image description here

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 *)

enter image description here

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 *)

enter image description here

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 *)

enter image description here

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 *)

enter image description here

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 *)

enter image description here

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 *)

enter image description here

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 *)

enter image description here

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

share|improve this answer
    
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 at 16:08
    
+1 extremely instructive and fun reading. How did you discover the EllipticE-alternative? Hopefully not with MMA :) –  eldo Aug 12 at 20:02
1  
@ 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 at 6:50
    
I am pleased you chose this site as your new "home" for such insights. Thank you! –  Mr.Wizard Aug 15 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 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 *)

enter image description here

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 *)

enter image description here

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 *)

enter image description here

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 *)

enter image description here

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.

share|improve this answer
    
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 at 15:16
    
@ Raffaele: please present your solution here so we can discuss it. Thank you. –  Dr. Wolfgang Hintze Aug 20 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

share|improve this answer
    
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 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]);
share|improve this answer
    
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 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 at 10:42
    
Sorry Wolfang, check the email. –  Raffaele Carlone Aug 28 at 11:01

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.