Trying to write a code to approximate a certain singular integral.

9 次查看（过去 30 天）
Lorenzo Wagemaker
Lorenzo Wagemaker 2016-9-28
编辑： John D'Errico 2016-9-29
Hey guys,
I'm fairly new at Matlab but would like to use it to compute the integral from [0,infinity): (sin(x^5)) / (x^(2)*(1+x)^(55)) dx
I realise that this can easily be solved with: fun = @(x) (sin(x.^5)) ./ (x.^(2).*(1+x).^55);
q = integral(fun,0,inf)
But would like to basically create a code that solves this. I was initially planning on using the Gaussian quadrature but quickly realised that it only works for polynomials. Not sure where to go from here, anyone know any other techniques/methods to solve such singular integrals which can be put in Matlab?
Thanks a lot
  3 个评论
显示 1更早的评论
Walter Roberson
Walter Roberson 2016-9-28
https://en.wikipedia.org/wiki/Gaussian_quadrature
"Gaussian quadrature as above will only produce good results if the function f(x) is well approximated by a polynomial function within the range [−1, 1]. The method is not, for example, suitable for functions with singularities. However, if the integrated function can be written as f ( x ) = ω(x)g(x) where g(x) is approximately polynomial and ω(x) is known, then alternative weights wiand points xi that depend on the weighting function ω(x) may give better results, where [...]"
John D'Errico
John D'Errico 2016-9-29
But that does not mean that Gaussian quadrature is only for polynomials.
For example, many (if not most) numerical quadrature methods will fail for problems with singularities. And it is also true that Gaussian quadrature is often a very good method for problems with specific known singularity types, since the singularity can then be resolved using an appropriate weight function.

请先登录，再进行评论。

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

回答（1 个）

John D'Errico
John D'Errico 2016-9-28
Whats the problem? Did you try it?
syms x
fun = (sin(x.^5)) ./ (x.^(2).*(1+x).^55)
I = int(fun,[0,inf]);
vpa(I)
ans =
0.00000079051133068352607021349297955875
Or using purely numerical methods, we can stop around 1 or so,because the large powers of x in that fraction ensure it underflows for x about 1 or 2.
fun = @(x) (sin(x.^5)) ./ (x.^(2).*(1+x).^55);
fun([.5 1 2 3])
ans =
2.5812e-11 2.3356e-17 7.9024e-28 -7.6182e-35
So an upper limit of 2 seems entirely far enough to get full precision.
format long g
integral(fun,0,2)
ans =
7.90511330683526e-07
  2 个评论
Lorenzo Wagemaker
Lorenzo Wagemaker 2016-9-28
Thank you for your reply, This is helpful, though my essay is mainly about the techniques involved in calculating these types of integrals. Thus simply compiling this short amount of code won't really do much for me. Instead I need to demonstrate what method matlab uses to calculate such an approximation. Do you happen to know hat quadrature rule matlab uses for such a computation?
Thank you
John D'Errico
John D'Errico 2016-9-29
编辑：John D'Errico 2016-9-29
That was not your question. Why do people ask one question, but really have something completely different in mind?
The symbolic computation I did computed the integral symbolically (nasty looking at that), then converted it to a numeric value.
The integral function, if you look at the documentation, has a reference.
[1] L.F. Shampine "Vectorized Adaptive Quadrature in MATLAB®," Journal of Computational and Applied Mathematics, 211, 2008, pp.131–140.
To be honest, I don't think that won't really help you that much, since what matters on these things is understanding the methods of numerical analysis that are available. There are many tools one should have available in your computational arsenal. A thorough understanding of them all, including Gaussian quadrature, adaptive methods, symbolic methods, etc., will give you the ability to solve problems. But no one method will suffice for all problems, and any specific method can often be tripped up by a careful choice of problem.

请先登录，再进行评论。

类别

在 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