Thread Subject: Integration Trouble

I'm trying plot the Kac Formula and I am having trouble. Here is a link to the formula, http://mathworld.wolfram.com/KacFormula.html (the formula that goes from 0 to 1).
It starts with an integral and I can't get my code to work. I created an m-file with this

syms t n

Integ=sqrt(1/(t^2-1)^2-((n+1)^2*t^(2*n))/(1-t^(2*n+2))^2);

int(Integ,t,0,1)

Matlab is just returning the integrand. Where is my problem?

Alan wrote:
> I'm trying plot the Kac Formula and I am having trouble. Here is a link
> to the formula, http://mathworld.wolfram.com/KacFormula.html (the
> formula that goes from 0 to 1).
> It starts with an integral and I can't get my code to work. I created
> an m-file with this
>
> syms t n
>
> Integ=sqrt(1/(t^2-1)^2-((n+1)^2*t^(2*n))/(1-t^(2*n+2))^2);
>
> int(Integ,t,0,1)
>
> Matlab is just returning the integrand. Where is my problem?


Well, quite simply, not every integral that can be framed has a
closed-form solution.

Walter Roberson wrote:
> Alan wrote:
>> I'm trying plot the Kac Formula and I am having trouble. Here is a
>> link to the formula, http://mathworld.wolfram.com/KacFormula.html (the
>> formula that goes from 0 to 1).
>> It starts with an integral and I can't get my code to work. I
>> created an m-file with this
>>
>> syms t n
>>
>> Integ=sqrt(1/(t^2-1)^2-((n+1)^2*t^(2*n))/(1-t^(2*n+2))^2);
>>
>> int(Integ,t,0,1)
>>
>> Matlab is just returning the integrand. Where is my problem?

> Well, quite simply, not every integral that can be framed has a
> closed-form solution.

If you examine Integ near t=1, you will see that at t=1 exactly, the
expression has two places in which there is a division by 0, once inside the
sqrt(), and the second place in the denominator. The first one, inside the
square root, factors as 1/((t-1)*(t+1))^2; because t is <= 0, the t-1 part is
<= 0 making t^2-1 <= 0, but the ^2 afterwards makes the expression >= 0. None
the less, we can see that this is not a removable discontinuity at this level:
there is no balancing expression in t in the numerator, so as t approaches 1,
the first term of the sum inside the sqrt() is going to grow to infinity. You
have to examine the integrand as a whole to determine what is happening near
1, and that in turn is going to depend upon the value of n -- in particular,
if n is not a positive integral multiple of 1/2 then you run into complexities
analyzing the powers, especially if n is complex -- and how is int() to know
that n is *not* complex?

With the information about n underdetermined, and with the singularity as
t->1, it is not surprising that a generalized integrator has difficulty
finding a closed form solution for this.


If you try work-arounds such as taylor() to get an approximation of the upper
and lower bounds of Integ (before you integrate), then although taylor()'ing
near t=1 looks complicated, it turns out in practice to evaluate to a simple
expression. However, if you try to taylor near 0, you will find that the first
differential of Integ is discontinuous at 0, so no taylor series can be formed
around t=0 . I was working on taylor'ing near 1/1000 but it was taking too
long even for a very modest number of terms.

In short... the expression is not readily integrated. There might be tricks to
integrate it... I'm not a calculus specialist.

Interestingly, the agreement with your reference for the first n 6 values is very good using a numerical approach.



Integ = @(t,n) sqrt(1./(t.^2-1).^2-((n+1).^2*t.^(2*n))./(1-t.^(2*n+2)).^2);
E = zeros(1,6);

for ii = 1:6
    E(ii) = quadgk(@(x) Integ(x,ii),0,.9999985,'abstol',1e-10)*4/pi;
end

x = [1 1 2 2 3 3 4 4 5 5 6 6 7];
y = reshape(repmat(E,2,1),1,[]);
line(x,[0 y])
xlim([0,7])
ylim([0,2])