Numerical integration of a function with a singularity

improper-integralsintegrationnumerical methodsnumerical-calculus

I'm trying to compute numerically a function like the following:
$$
F(t)=\int_{0}^{t}{\frac{f(\tau)}{\sqrt{t-\tau}}d\tau}
$$
I tried to adopt the composite Simpson's rule, but the problem is that when $\tau=t$, I have a division by $0$, so I can't evaluate the integrand at the last point. I have also tried some substitutions but without success. Could someone please help me with how can I proceed to get rid of it? Thank you in advance!

Best Answer

If you consider the change of variable $t-\tau = y^2$ the integral becomes $$ \int_{\sqrt t}^0 \frac{f(t-y^2)}{y}\cdot (-2y) dy = 2 \int_0^{\sqrt{t}} f(t-y^2) dy. $$

Assuming that the singularity was expressed in the term $\frac{1}{\sqrt{t-\tau}}$ and $f$ is nice enough, Simpson's method will now work in standard way.

Related Solutions

[Math] Numerical integration in Matlab (Simpson’s rule)

A couple of notes. Since the integrand is never less than $1$, you know that $x_l<170$. In fact, when $x=170$, the integrand will be about $1.087$, so the eventual answer is bounded below by $\frac{170}{1.087}=156.5$. That gives you a pretty good idea about the starting point.

Bisection is useful when you don't have the derivative available, but here the you have the derivative in the form of the integrand, so it seems to be simpler to implement Newton's method in this case unless bisection is the point of the assignment.

Also you know the integral up to a point near to the actual root at each step after the first, so you need not integrate all the way from zero to the next guess for the root, just from the previous guess. Even if you integrated from zero every time, this program would zip through pretty fast anyway.

EDIT: Newton's method

% simp.m

f = @(x) sqrt(1+(x^2/68000)^2);

err = 1;
tol = 1.0e-8;
N = 10000;
a = 0;
b = 163;
while abs(err) > tol,
    h = (b-a)/N;
    y = f(a);
    for i = 1:N/2,
        y = y+4*f((2*i-1)*h)+2*f(2*i*h);
    end
    y = y+4*f(b-h)+f(b);
    y = h/3*y-170;
    yp = f(b);
    err = y/yp;
    b = b-err;
end
b

[Math] Evaluating integral with a singularity.

You have two problems.

The problem you ask about, why doesn't this code work, belongs somewhere other than MSE; it's not a math question. (When you ask this question again elsewhere you should say what the error message is...)

There's also confusion with the math. Strictly speaking there's no such thing as $$\int_{-1}^1\frac{dt}{t};$$you need to decide how to "interpret" this intergral. The most common interpretation would be $$\lim_{\epsilon\to0^+}\int_{\epsilon<|t|<1}\frac{dt}t,$$known as the "principal value" integral. This is not the same as $$\lim_{\epsilon\to0^+}\int_{-1}^1\frac{dt}{t+i\epsilon}.$$In fact here the principal-value integral is $0$, while that complex trick gives $i\pi$.

Which one is right? Neither. Because neither is actually equal to $\int_{-1}^1dt/t$; that integral simply does not exist. Which one do you want? That depends on why you think you want to find this integral in the first place. (But when one writes $\int_{-1}^1dt/t$ without expplanation, the most common interpretation is that the author meant the principal value, which again is not what that trick you read in a book gives.)

Best Answer

Related Solutions

[Math] Numerical integration in Matlab (Simpson’s rule)

[Math] Evaluating integral with a singularity.

Related Question