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

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

[John D'Errico]

I had to laugh at the claim that Gaussian quadrature ONLY works for polynomials. Did you read that somewhere? Then what you read was flat out wrong. Did you somehow draw that conclusion from something that you read? Again, wrong.

Is Gaussian quadrature appropriate for this problem? Probably not. At the same time, don't do yourself and anyone else who reads that claim a disservice.

[Walter Roberson]

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]

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.

Answer 1

[John D'Errico]

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

[Lorenzo Wagemaker]

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]

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.