Numerical integration with a singularity at the upper limit

How can I numerically integrate a function with a singularity at the upper limit? The 'integral' function is not smooth to changes in the lower limit. Here is a code with my function and a plot that illustrates the instability:

n=20;
h1=@(x) x.^(1/2+n).*(1-x).^(1/2-n);
inth1=@(p) integral(h1,p,0.99);
fplot(inth1,[0.5,0.55])

quadgk does not solve the problem.

Answer 1

[John D'Errico]

No matter what numerical methods scheme one chooses to solve this problem, you can always choose a sufficiently large value of n that makes it unstable. NO numerical method will survive your test, IF you are willing to make your function arbitrarily nasty.

Anyway, there is not in fact a singularity at the limits of integration here, but beyond those limits. The singularities are at x=0 and x=1, whereas you have chosen [p,0.99] as the limits of integration. This is a completely different problem than one with a singularity.

If you insist on solving such arbitrarily nearly singular problems, you will probably need to work in a higher precision arithmetic than double precision. So I'd start by looking at the symbolic toolbox.

[Linas Tarasonis]

Thank you, John, for your answer. The symbolic math toolbox works well if n is not too large indeed.

Answer 2

[Torsten]

Your integral only exists for n < 3/2 (if the lower x-limit is greater 0).

Best wishes

Torsten.

[Linas Tarasonis]

Thank you for your answer, Torsten. This is why the upper limit is 0.99 instead of 1.

[Torsten]

0.99 is somewhat arbitrary if the integral does not exist, isn't it ?

Best wishes

Torsten.