I cannot make a calculation with the command int

Hello, I am doing a calculation, but matlab doesn't want to do it with the command int. I have found the solution with other software or using the command trapz, but it is linger and less accurate.

this is the code

syms x

a=int((72290819601*x^2*(x - 1)^10*(x + 1)^10*(261*x^4 - 90*x^2 + 5)^2)/(134217728*(x - 4752731968273649/18014398509481984)^2),x,-1,1)

It gives infinity as result, whereas it is equal to 3831.098 (all the digits in my formula comes from a simplify command in matlab.

I am running Matlab 2010b

Thanks a lot for any suggestion for fixing the problem, otherwise I'll use the longer procedure from trapz

[Sean de Wolski]

+1, good question in both content and description!

Answer 1

[Honglei Chen]

looks to me you have a singular point between -1 and 1 at 4752731968273649/18014398509481984. This might be why

[mortain]

It's not really a singular point...I specify at the third answer

[John D'Errico]

How is it not singular? The divide by zero works for you? Wow. It must be nice to live in a world where 1/0 has a finite value. I wonder when that changed? The numerator is not zero at that point. Sorry. It IS singular.

[mortain]

I said it, because if you plot the graph, it does not seem to have a discontinuity point there. In addition, the numerator, according to what a friend told me, is zero as well at that point. However, I'll check again. I would love a world where1/0 has finite value, but then it would be a mess :-) !

[Walter Roberson]

You are not plotting at a fine-enough resolution to see the infinity.

[mortain]

I tried to see really near the point...but I could not find any asymptote...you can give a try on any online software or matlab itself....I reckon the formula I put carry a certain round off error, hence the numerator is not really equal to zero, but E-20, which is comparable to a round off error. I would like to have an analytical solution of the integral. Thanks for the suggestions.

Can you please tell me how I can see the infinity or more digits, please, I tried using the digits(40) at the beginning, but I could not.

Answer 2

[Walter Roberson]

Try increasing your Digits setting. The integral has some large intermediate terms.

[mortain]

What do you mean by increasing the digits setting?

[Walter Roberson]

http://www.mathworks.com/help/toolbox/symbolic/digits.html

[mortain]

Tomorrow in the lab I'll try. So you suggest me to write at the beginning digits(40) at the beginning.

However, as I mentioned in the third reply, matlab assumes that that point (0.26 ish) gives infinite which is not true.

[Walter Roberson]

Sorry, I had a disk problem at home a few days ago and had to wipe my disk; I haven't reinstalled my software yet, so I cannot test at the moment.

Answer 3

[Sean de Wolski]

When x = 4752731968273649/18014398509481984 or 0.26ish, the denominator is zero, the numerator is not zero and the value is infinite.

In the numerical world, it would fail if this point was evaluated. You just don't have this in your trapz calculation:

Let's make the function in the happy numerical world:

f = str2func(vectorize('@(x)(72290819601*x^2*(x - 1)^10*(x + 1)^10*(261*x^4 - 90*x^2 + 5)^2)/(134217728*(x - 4752731968273649/18014398509481984)^2)'));

And then evaluate it at our trouble point:

f(4752731968273649/18014398509481984)

Boom!

[mortain]

Hello Sean,

Thanks for the like. However, the point is actually not a boom point. A friend showed that if you apply Hopital couple of times, you get a value for that point. It is actually a 0/0 and apparently matlab is not able to cope with it.

I think I'll stick with trapz or with a sum, even though I would have preferred to get a great accuracy by using not the numerical analysis.

Thanks for the answer and the time, I guess there's no way for gettin an analytical solution.

[John D'Errico]

Just because trapz gives a finite result, this does not make that result correct. You should recognize that this integral has no finite result. int is telling you that. Listen to what it says.

[John D'Errico]

Mortain, your friend is wrong here. Since the numerator is not truly zero at that point, it is not 0/0. A really small value does not count. At that location, the numerator evaluates to (roughly) 0.000000000000000000063950710191408455365893342257719, not zero.

[Walter Roberson]

http://www.wolframalpha.com/input/?i=Limit[%2872290819601*x^2*%28x+-+1%29^10*%28x+%2B+1%29^10*%28261*x^4+-+90*x^2+%2B+5%29^2%29%2F%28134217728*%28x+-+4752731968273649%2F18014398509481984%29^2%29%2Cx-%3E4752731968273649%2F18014398509481984]

The result there is infinity.

[Richard Brown]

I can add Maple to the list of software packages confirming this :)

[mortain]

Sure guys that it is not the fact that I have only some digits which brings a certain value of the numerator, so it is more related to a round off error.

I add the code here at the bottom.

Answer 4

[mortain]

    close all
    clear all
    clc
    n=5
    alfa=10;
    beta=10;
    digits(40)
    syms x
dim=5; %max order polynomial
b=zeros(dim,dim); %initialization of the roots matrix
df=simplify((-1)^n/(2^n*factorial(n))*diff((1-x)^(alfa+n)*(1+x)^(beta+n),n)/((1-x)^(alfa)*(1+x)^(beta)));

%creation of the orthogonal polynomial of order n
df=simplify((-1)^n/(2^n*factorial(n))*diff((1-x)^(alfa+n)*(1+x)^(beta+n),n)/((1-x)^(alfa)*(1+x)^(beta)));
%Creation of the coefficients for the matrix
a=sym2poly(simplify((-1)^n/(2^n*factorial(n))*diff((1-x)^(alfa+n)*(1+x)^(beta+n),n)/((1-x)^(alfa)*(1+x)^(beta))));
%dummy vector for storing all the roots in one matrix
asd=zeros(1,dim);
bbb=roots(a);
for l=1:n
    asd(l)=bbb(l);
end
b(n,:)=asd;
for i=1:n
%Evaluation of the first derivative of the polynomial expression
dff=diff(df,1);
%Evaluation of the first moltiplicando of the weight calcultion
df3(n,i) =1/(subs(dff, 'x', b(n,i)))^2;

%Evaluation of the integrand fuction=integ
integ(n,i)=(gamma(alfa+beta+2)*(1-x)^alfa*(1+x)^beta/(2^(alfa+beta+1)*gamma(alfa+1)*gamma(beta+1)))*((df)/(x-b(n,i)))^2

% Numerical evaluation of the integral
integr(n,i)=int(integ(n,i),x,-1,1)
end

In addition I am doing some other integration and Matlab is not able to integrate this

int((7186705221432913*(x^2 - 1)^2)/(18014398509481984*exp(x^2/2)*(x + 1)^2), x = -Inf..Inf)/4

this time saying that Warning: Explicit integral could not be found.

I am actually looking for an explicit formulation, but if it is impossible a numerical approximation. I guess the trapezoidal rule is embedded already, or you have other suggestions?

Thanks

[Walter Roberson]

For that last integral, Maple indicates the answer is

(7186705221432913/36028797018963968)*sqrt(Pi)*sqrt(2)

[bym]

FWIW, in matlab (I just learned the simplify 'trick' for integrals)

int((7186705221432913*(x^2 - 1)^2)/(18014398509481984*exp(x^2/2)*(x + 1)^2), -Inf,Inf)/4

Warning: Explicit integral could not be found.

> In sym.int at 64

ans =

int((7186705221432913*(x^2 - 1)^2)/(18014398509481984*exp(x^2/2)*(x + 1)^2), x = -Inf..Inf)/4

>> simplify(ans)

ans =

(7186705221432913*2^(1/2)*pi^(1/2))/36028797018963968

Answer 5

[Walter Roberson]

> Digits := 50;
                               50
> plot(72290819601*x^2*(x-1)^10*(x+1)^10*(261*x^4-90*x^2+5)^2/(134217728*(x-4752731968273649/18014398509481984)^2), x = .2638296 .. .2638297);

The singularity is clearly present, and with 50 digits of precision being calculated, it sure isn't just round-off error in calculating the plot.

Perhaps you are trying to say that round-off error in your generating calculations led to an incorrect formula that you then asked us to integrate? If so then what is the correct formula?