Got Questions? Get Answers.
Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Integration Trouble

Subject: Integration Trouble

From: Alan

Date: 16 Mar, 2010 03:07:05

Message: 1 of 4

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?

Subject: Integration Trouble

From: Walter Roberson

Date: 16 Mar, 2010 05:49:54

Message: 2 of 4

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.

Subject: Integration Trouble

From: Walter Roberson

Date: 16 Mar, 2010 23:53:31

Message: 3 of 4

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.

Subject: Integration Trouble

From: Matt Fig

Date: 17 Mar, 2010 03:57:07

Message: 4 of 4

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

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us