I cannot make a calculation with the command int

5 次查看（过去 30 天）
mortain
mortain 2012-5-2
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

请先登录，再回答此问题。

回答（5 个）

Honglei Chen
Honglei Chen 2012-5-2
looks to me you have a singular point between -1 and 1 at 4752731968273649/18014398509481984. This might be why
  5 个评论
显示 3更早的评论
Walter Roberson
Walter Roberson 2012-5-3
You are not plotting at a fine-enough resolution to see the infinity.
mortain
mortain 2012-5-3
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.

请先登录，再进行评论。

Walter Roberson
Walter Roberson 2012-5-2
Try increasing your Digits setting. The integral has some large intermediate terms.
  4 个评论
显示 2更早的评论
mortain
mortain 2012-5-3
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
Walter Roberson 2012-5-3
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.

请先登录，再进行评论。

Sean de Wolski
Sean de Wolski 2012-5-2
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!
  6 个评论
显示 4更早的评论
Richard Brown
Richard Brown 2012-5-3
I can add Maple to the list of software packages confirming this :)
mortain
mortain 2012-5-3
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.

请先登录，再进行评论。

mortain
mortain 2012-5-3
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
  2 个评论
Walter Roberson
Walter Roberson 2012-5-11
For that last integral, Maple indicates the answer is
(7186705221432913/36028797018963968)*sqrt(Pi)*sqrt(2)
bym
bym 2012-5-11
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

请先登录，再进行评论。

Walter Roberson
Walter Roberson 2012-5-11
> 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?

类别

在 Help Center 和 File Exchange 中查找有关 Numerical Integration and Differentiation 的更多信息

标签

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by