I'm looking for $I=\int_{0}^\infty f(u)du$ where $f$ is the solution to the integral equation $$\frac{\pi}{2}f(u) = \frac{1}{1+u^2} - \int_{0}^\infty \frac{f(v) dv}{4+(u-v)^2}$$ for $u$ between $0$ and $\infty$.

While knowing $f(u)$ exactly would be nice, I'd be quite happy with just $I=\int_{0}^\infty f(u)du$ or an approximation to $I$! My goal would be an approximation within $1$%.

I have two attempts below - one with an ansatz for $f(u)$ and the other a perturbative series for $I$. I'm hopeful that either direction will be fruitful.

Here's a sketch of my attempt that roughly estimates $I \approx .6$.

First, note that a problem with different bounds and with $u$ between $\infty$ and $\infty$, $$\frac{\pi}{2}g(u) = \frac{1}{1+u^2} - \int_{\color{red}{-\infty}}^\infty \frac{g(v) dv}{4+(u-v)^2}$$ can be solved via Fourier transform to find $g(u) = \frac{1}{e^{\frac{\pi}{2}u}+e^{-\frac{\pi}{2}u}}$ and $\int_{0}^\infty g(u) du = \frac{1}{2}$.

I anticipate that $I$ for the problem of interest will be a little more than $.5$, given that it appears we're subtracting less from the right hand side when we use smaller bounds.

Consider the ansatz $$f(u) = \frac{c_1}{e^{\frac{c_2 \pi}{2}u}+e^{-\frac{c_2\pi}{2}u}}$$ Changing the parameters $c_1$ and $c_2$ by hand, I find with numerical integration that $c_1=1.11$ and $c_2=.93$ make the left hand side and right hand side of $$\frac{\pi}{2}f(u) = \frac{1}{1+u^2} - \int_{0}^\infty \frac{f(v) dv}{4+(u-v)^2}$$ nearly equal as a function of $u$. I anticipate this ansatz won't include the exact solution, but it gives what I think is a reasonable estimate for $I$ of about $.6$.

Here is another attempt. Note that these linear problems can in principle be solved through perturbative solutions; $$f(u) = h(u)+\lambda \int_0^\infty K(u,v) f(v)$$ is solved by $$f(u) = \sum_{n=0}^\infty \lambda^n \int_0^\infty du_2 ... \int_0^\infty du_{n+1} K(u, u_2)...K(u_n, u_{n+1}) h(u_{n+1})$$ so $$I = \sum_{n=0}^\infty \lambda^n \int_0^\infty du_1 \int_0^\infty du_2 ... \int_0^\infty du_{n+1} K(u_1, u_2)...K(u_n, u_{n+1}) h(u_{n+1}).$$

That is, we have $$I = \sum_{n=0}^\infty \left(-\frac{2}{\pi}\right)^n \frac{2}{\pi} \int_0^\infty du_1 ... \int_0^\infty du_{n+1} \frac{1}{4+(u_1-u_2)^2}...\frac{1}{4+(u_n-u_{n+1})^2}\frac{1}{1+u_{n+1}^2}$$

The first few partial sums are $1, 0.318, 0.866$. Assuming the terms are decreasing in the series, the alternating sum is then bounded between $0.318<I<0.866$. This is consistent with $I\approx.6$, but the bounds are quite large, and the integrals rapidly become both analytically and computationally challenging to evaluate.

  • $\begingroup$ And the motivation is? Is this an exercise known to be solvable, a research task, a curious question you dreamt up one day... ?# $\endgroup$
    – FShrike
    Commented Mar 17, 2023 at 18:05
  • $\begingroup$ @FShrike These types of integral equations appear in describing the ground states of Bethe ansatz-solvable Hamiltonians; the version with bounds going down to $-\infty$ describes the occupation numbers of the spin-1/2 Heisenberg chain in the thermodynamic limit and be used to derive the Bethe-Hulthen energy value. $\endgroup$
    – user196574
    Commented Mar 18, 2023 at 23:22
  • 2
    $\begingroup$ Your equation is a particular case of the generalised Love equation. Numerical approaches are discussed for instance in Farina et al. (2020). $\endgroup$
    – TheSimpliFire
    Commented Mar 19, 2023 at 11:36
  • 1
    $\begingroup$ @TheSimpliFire I looked at Farina's paper. There is no hint there whatsoever of how to deal with slowly decaying functions. Also the "numerical" section lacks any details about how to estimate the approximation errors rigorously for any of the proposed methods. Utter disappointment! $\endgroup$
    – fedja
    Commented Mar 20, 2023 at 16:08
  • 1
    $\begingroup$ @fedja Oh, that's a shame. The good thing is, now that I know the terminology I should be able to find more useful references. $\endgroup$
    – TheSimpliFire
    Commented Mar 20, 2023 at 18:38

2 Answers 2


This question would be better suited for MO, IMHO.

Here is a quick and dirty way to do it at the engineering level of rigor. There are two questions here: one is how to compute the integral with decent precision and another one is how to justify that we are not talking nonsense from the very beginning (It is obvious that $f\in L^2$ but there seems to be no easy way to show that it is in $L^1$).

Let's write the equation as $f+Tf=\varphi$ where $(Tf)(u)=\int_0^\infty K(u-v)f(v)\,dv$, $K(t)=\frac 2\pi\frac 1{4+t^2}$ (so $\int_\mathbb R K=1$), and $\varphi(t)=\frac 2\pi \frac 1{1+t^2}$.

The first observation is that $T$ is self-adjoint ($K$ is even) positive definite $\widehat K>0$ operator of norm $\le 1$ ($\|K\|_{L^1}=1$) in $L^2=L^2(0,+\infty)$. Thus its spectrum lies on $[0,1]$, so $I+T$ is, indeed, invertible, and $\|(I+T)^{-1}\|_{L^2\to L^2}\le 1$.

We shall now solve the auxiliary equation $g+Tg=1$. That seems to make little sense from the $L^2$ point of view, but if we write $g=\frac 12+h$, then it is equivalent to $h+Th=\frac 12(1-T1)=\psi$ and $\psi$ is now in $L^2$ (it is $\psi(u)=\frac 12\int_u^{\infty}K(v)\,dv$, so it decays like $1/u$ at infinity.

Thus, we can write $$ \langle f,1\rangle=\langle f, g+Tg\rangle= \langle f+Tf,g \rangle=\langle \varphi, g\rangle \\ = \frac 12\langle\varphi,1\rangle+\langle \varphi,h\rangle =\frac 12+\langle \varphi,h\rangle $$ and since $\varphi\in L^2$, it suffices to know $h$ with high precision in $L^2$ to find $\langle f, 1\rangle$ with high precision regardless of how large $\|f\|_{L^1}$ is as long as it is finite.

Now we are just solving $(I+T)h=\psi$. We want to approximate $(I+T)^{-1}$ by a polynomial $P(T)$ of low degree. The $L^2\to L^2$ norm of the difference will be just $E(P)=\max_{t\in[0,1]}\left|P(t)-\frac 1{1+t}\right|$. What you tried to do was to take $P(t)=1-t+t^2-\dots\pm t^n$ which makes as much sense as trying to compute $\pi$ directly from the Leibniz formula with the alternating harmonic series. You can play a bit yourself with finding the best polynomial approximation, but just using the interpolation polynomial of degree $4$ with Chebyshev nodes gives you the trivial guaranteed bound of $\frac 1{256}$. Thus, using that approximation, you'll be off by at most $\frac 1{256}\|\psi\|_{L^2}\|\varphi\|_{L^2}$.

Since $\|\varphi\|_{L^2}=\frac 2\pi\sqrt{\frac\pi 4}=\frac{1}{\sqrt{\pi}}$ and $$ \psi(u)=\frac 1\pi\int_u^\infty\frac{dv}{4+v^2}\le\frac 1\pi\int_u^{\infty}\frac{2\,dv}{(2+v)^2}=\frac 2\pi\frac{1}{2+u}\, $$ we get $\|\psi\|_{L^2}^2\le \frac 2{\pi^2}$, so $\|\psi\|_{L^2}\|\varphi\|_{L^2}\le \sqrt{\frac 2{\pi^3}}$ yielding the garanteed upper bound for the error of $0.001$ plus all the errors you accumulate when computing $T$ and $T^2$ of various functions numerically (there is no real need to try to compute $T^4$ because $\langle T^4 \alpha, \beta\rangle=\langle T^2\alpha, T^2\beta\rangle$.

Now comes the second issue: the functions involved are all decaying slowly, so to find a decent numeric integration scheme is not so simple. I hope Farina's book has something useful to say about that, but if not, then we'll have to discuss this point too, though at the moment I do not have any much brighter idea than using Simpson on a reasonably fine lattice stretching far enough. I'll try some numerics to see what makes sense and what doesn't but it will take some time :-)

Hope this helps a bit.

Edit Here is the missing piece: the proof that the integral, indeed, exists. We shall show that $I+T$ is invertible in $L^2$ with the weight $\alpha(u)=\beta(u)^2$ where $\beta(u)=1+cu^p$ with any $\frac 12<p< 1$ and sufficiently small positive $c<c(p)$. To this end, we just need to show that the self-adjoint part $S$ of $T$ in that space satisfies $\delta I+S\ge 0$ with some $\delta<1$. It is equivalent to showing that the self-adjoint operator $$ Sf(u)=\int_0^\infty \widetilde K(u,v)f(v)\,dv $$ with $$ \widetilde K(u,v)=\frac 12\left[\frac{\beta(u)}{\beta(v)} \frac{\beta(v)}{\beta(u)}\right] K(u,v) $$ is almost positive definite in the usual non-weighted $L^2$. We shall just compare it with $K$ and write the difference as the integral operator with the kernel $$ \widetilde K(u,v)=\frac 12\left[\frac{\beta(u)}{\beta(v)}- \frac{\beta(v)^2}{\beta(u)^2}\right]^2 K(u,v) $$ All we want is to estimate the norm of this operator by a small constant when $c$ is small. Invoking the Schur test with the function $\psi(u)=\frac{1}{\beta(u)}$, we see that it suffices to show that $$ \frac 12\int_0^{\infty} \frac{\beta(u)}{\beta(v)}\frac 12\left[\frac{\beta(u)}{\beta(v)}+ \frac{\beta(v)}{\beta(u)}-2\right] K(u,v)\,dv $$ is small uniformly in $u>0$.

Note that $\frac 12(Z+Z^{-1}-2)\le\max(Z,Z^{-1})-1$, so we need just to estimate $$ I_1=\int_u^{\infty} \frac{\beta(u)}{\beta(v)}\frac 12\left[ \frac{\beta(v)}{\beta(u)}-1\right] K(u,v)\,dv $$ and $$ I_2=\int_0^{u} \frac{\beta(u)}{\beta(v)}\frac 12\left[ \frac{\beta(u)}{\beta(v)}-1\right] K(u,v)\,dv\,. $$

The estimate of $I_1$ is easy: just note that $\frac{\beta(x)}{\beta(x+y)}$ is an increasing function of $x$ for fixed $y>0$, so $$ I_1\le\int_u^\infty\left(1-\frac{\beta(0)}{\beta(v-u)}\right)K(u,v)\,dv =\frac 2\pi\int_0^\infty\left(1-\frac 1{1+ct^p}\right)\frac{1}{4+t^2}\,dt $$ and the dominated convergence theorem implies that this bound tends to $0$ as $c\to 0$.

$I_2$ is trickier. We'll split it further into $I_2'+I_2''=\int_0^{u/2}+\int_{u/2}^u$.

For $I_2''$, we have $\frac{\beta(u)}{\beta(v)}\le 2^p$, so $$ I_2'\le 4^p\int_{u/2}^u\left(1-\frac{\beta(v)}{\beta(u)}\right)K(u,v)\,dv\le 4^p\int_{u/2}^u\left(1-\frac{\beta(0)}{\beta(u-v)}\right)K(u,v)\,dv $$ and we finish just as for $I_1$.

For $I_2'$, we have $K(u,v)\le \frac 1{4+\frac{u^2}4}$. Now we first notice that $\beta(u)/\beta(v)$ is increasing with $c$. However, for $c=1$, we have $$ I_2'\le\frac{(1+u^p)^2}{{4+\frac{u^2}4}}\int_0^\infty \frac{dv}{(1+v^p)^2} $$ and since the integral converges and $2p<2$, we conclude that $I_2'$ is small for sufficiently large $u>u_0$ where $u_0$ does not depend on $c$. But for $0<v<u<u_0$ the integrand tends to $0$ uniformly as $c\to 0$, so we are done.

The challenge now becomes to figure out if $f$ is, indeed, non-negative.

  • $\begingroup$ @TheSimpliFire Actually, the quadratic polynomial $P(t)=0.9934-0.8296t+0.3428t^2$ approximates $\frac 1{1+t}$ with precision $0.008$ already, which, multiplied by $\sqrt{\frac 2{\pi^3}}$ is less than $0.0021$, which is less than half-percent of $1/2$. Also $T\varphi$ and $\psi$ are elementary functions. So, only $\langle T^2\psi,\varphi\rangle$ is really hard to find with high precision. I'll think of how that can be done now. $\endgroup$
    – fedja
    Commented Mar 20, 2023 at 15:38
  • $\begingroup$ @TheSimpliFire Since $\psi(u)=\frac 1\pi(\frac \pi 2-\arctan \frac u2)$, we can turn it into a rational function by using the approximation $\frac\pi 2-\arctan x\approx\frac{\frac\pi 2+1.75x+x^2}{(1+x)(1+0.75x+x^2)}$ with uniform error $<0.0016$, so the uniform error in $T\psi$ is bounded by $1/\pi$ times that and we also have $\|\varphi\|_{L^1}=1$, so $\|T\varphi\|_{L^1}\le 1$ as well and the error in $0.3428\langle T\psi,T\varphi\rangle$ will be $<0.0002$, so we have already used $0.0023$ and with the OP request for $0.005$, we have $0.0027$ left on numeric integration of explicit functions. $\endgroup$
    – fedja
    Commented Mar 20, 2023 at 20:13
  • $\begingroup$ @TheSimpliFire With all that junk, an online symbolic integrator that produced ugly but precise antiderivatives of the rational functions involved, and brute force Simpson for explicit functions, I got $I=\int_0^\infty f(x)\,dx=0.616\pm 0.005$, which was sort of expected from the beginning. Probably the error is actually at least twice smaller, so I'm pretty sure that $0.61<I<0.62$, if the integral exists at all, which I still have no idea how to prove... Now let's challenge someone to get a few more digits :lol: $\endgroup$
    – fedja
    Commented Mar 20, 2023 at 23:27
  • 1
    $\begingroup$ @user196574 Should I wait until the bounty expires to mark this answer as accepted? If you ask me, I would rather wait for somebody to prove that $f\in L^1$. But the choice is yours, there are no strict or unwritten rules here :-) $\endgroup$
    – fedja
    Commented Mar 21, 2023 at 15:36
  • 1
    $\begingroup$ @Simplifire. OK, I figured out why $f\in L^1$. There seems to be no good way to do it directly, but you can show that it is in $L^2$ with the weight $(1+cu^p)^2$ for any $p<2/3$ and choosing $p$ above $1/2$, we conclude it by Cauchy-Schwarz. So, all estimates make sense, indeed :-). However, I still do not know whether $f$ is positive or not. $\endgroup$
    – fedja
    Commented Mar 24, 2023 at 2:38

Both of your attempts are interesting, and they provide some insight into the problem. However, they don't provide a guaranteed accurate approximation of $I$ within 1%. Here is a different approach, using numerical methods to solve the integral equation and obtain a value for $I$.

  1. Discretize the integral equation: We can approximate the integral equation by discretizing it and turning it into a system of linear equations. Divide the interval $(0, \infty)$ into $N$ subintervals $(u_1, u_2, ..., u_N)$, where $u_i = i \Delta u$ for some small $\Delta u$. Then, the integral equation can be approximated as: $$\frac{\pi}{2} f(u_i) \approx \frac{1}{1 + u_i^2} - \sum_{j=1}^N \frac{f(u_j) \Delta u}{4 + (u_i - u_j)^2} $$

  2. Now we have a system of $N$ linear equations in the unknowns $f(u_1), f(u_2), \dots, f(u_N)$. Let $A$ be the $N \times N$ matrix with elements $a_{ij}$, and let $x_i = f(u_i)$ and $b_i = \frac{\pi}{2}f(u_i)$. Then, we can rewrite the system of equations as $Ax = b$, where: $$a_{ij} = \begin{cases} \frac{\pi}{2} - \frac{\Delta u}{4+(u_i-u_j)^2}, & i \neq j \\ \frac{\Delta u}{4+(u_i-u_j)^2}, & i = j \end{cases}\\ b_i = \frac{1}{1+u_i^2} $$

To solve the linear system $Ax = b$, we can use numerical techniques such as the Gauss-Seidel method or the conjugate gradient method.

  1. Compute the integral: Once we have the values of $f(u_i)$, we can approximate the integral $I$ by a simple Riemann sum: $$I \approx \sum_{i=1}^N f(u_i) \Delta u$$ To obtain an approximation within 1%, you can experiment with different values of $N$ and $\Delta u$ until the desired accuracy is reached. Keep in mind that increasing $N$ will improve the accuracy of the approximation but also increase the computational cost.

To summarize, although an exact solution to the problem is not provided, the numerical approach allows you to obtain an approximation within the desired accuracy.

  • 1
    $\begingroup$ ":you can experiment with different values of N and Δu until the desired accuracy is reached" How exactly do you know when to stop to guarantee the desired accuracy and not just to start a new religious belief and that you'll have enough computational power to reach that stopping point? (because it is the guarantee that is difficult here, not the computation itself). Of course, your approach works "in principle" but I'm afraid that the pesky details with error estimates will be too many to make the justification of 1% accuracy convincing. $\endgroup$
    – fedja
    Commented Mar 25, 2023 at 19:35

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .