NDSolve Convergence Issue

Question

I am trying to use NDSolve to solve r[u] equation below. The problem seems like the slow convergence. I have searched some questions regarding this kind of problem. I think many people answered in some form of increasing working precision. If I increase the workingprecision, i get that the argument is not upto the machine precision. I will attach the warning it gave when I do NDSolve and my code.

NIntegrate::slwcon: Numerical integration converging too slowly; suspect one 
of the following: singularity, value of the integration is 0, highly 
oscillatory integrand, or WorkingPrecision too small.
NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 
18 recursive bisections in s near {t,s} = {4.85075,0.0000644015}. NIntegrate 
obtained 5.58828*10^-21-7.99344*10^-27 I and 2.9124298522935253`*^-24 for 
the integral and error estimates.

Code:

L[u_] := 1/(64 \[Pi]^2 t^5 Sqrt[u]) (-4 I t (s + I t^2 + s Log[4]) + (s^2 + 
4 I s t - 4 t^2 + t^4) Log[-2 + (I s)/t + t] - (s^2 - 4 I s t - 4 t^2 + t^4) 
Log[2 + (I s)/t + t])

H[u_?NumericQ] := NIntegrate[t^2*u^(5/2)/2*(Log[1 - L[u]] + L[u]/(1 - 
L[u])), {t, 0, 10^3}, {s, -10^3, 10^3}, AccuracyGoal -> 8, PrecisionGoal -> 
8,WorkingPrecision -> 12]

Plot[H[u], {u, 0, 100}, PlotRange -> All]

sol = NDSolve[{H[u] == u*r'[u] - 2 r[u], r[10^6] == -1}, r, {u, 
10^-9, 10}, AccuracyGoal -> 8, PrecisionGoal -> 8, WorkingPrecision -> 12]

Plot[Evaluate[r[u] /. sol], {u, 10^-9, 10}, PlotRange -> All, PlotStyle -> 
Green]

Where on the NDSolve line, I wanted a boundary condition that $r\rightarrow -1$ as $u\rightarrow\infty$, but just set u=10^6 to avoid infinity expression; I hope that's okay.

Any help will be greatly appreciated! Thank you so much for reading my question!

Answer 1

This started as an extended comment but has morphed into an answer.

H[u_?NumericQ] := NIntegrate[t^2*u^(5/2)/2*(Log[1 - L[u]] + L[u]/(1 - L[u])), 
    {t, 0, 10^3}, {s, -10^3, 10^3}]
Table[H[u]/u^3 /. u -> 10^i, {i, 6}]
(* {414.781, 56.4169, 7.66857, 0.893855, 0.141469, 0.0192026} *)

plus the first warning message given in the question. Increasing WorkingPrecision yields

H[u_?NumericQ] := NIntegrate[t^2*u^(5/2)/2*(Log[1 - L[u]] + L[u]/(1 - L[u])), 
    {t, 0, 10^3}, {s, -10^3, 10^3}, WorkingPrecision -> 30] // Chop
(* {414.780936186193432754748625707, 56.4169472731001172506090076002, 
    7.66856870136135201710072009702, 1.041808172243427818769297282494, 
    0.141469419050140253351266614811, 0.0192025581519389795391545651320} *)

plus both the first and second warning messages in the question. Nonetheless, the two sets of answers are in good agreement, suggesting that they are correct to several significant figures. Incidentally, Chop has deleted imaginary parts of order 10^-15. Although the integrand as written is complex, it is possible to rewrite the integral symbolically to make the integrand real, although the resulting expression is impractically large.

Next, the ODE can be solved symbolically in terms of H.

Clear[H]
Flatten@DSolve[H[u] == u*r'[u] - 2 r[u], r[u], u]
(* {r[u] -> u^2*C[1] + u^2*Integrate[H[K[1]]/K[1]^3, {K[1], 1, u}]} *)

Thus, for r[u] to approach -1 as u approaches infinity, the integral must approach C[1] - 1/u^2. However, from the sample evaluations of H[u] above, it is clear that H[u] itself approaches 0 more slowly than 1/u. Therefore, the integral blows up as u approaches infinity. So, the boundary condition sought in the question cannot be achieved.

Addendum - Behavior of H[u] near u == 0.

Experimentation suggests that the double integral defining H[u] behaves a bit better at small u with options, Exclusions -> {0, 0}, AccuracyGoal -> Infinity, PrecisionGoal -> 6, MaxRecursion -> 12. Then, H[u]/u^3 scales roughly as u^(-6/7}, as can be seen from

LogLogPlot[H[u]/u^3, {u, 10^-15, 10^-1}, PlotRange -> All, ImageSize -> Large, 
    AxesLabel -> {u, H}, LabelStyle -> Directive[12, Bold, Black]]

The integral can, then, be performed without error messages (as requested by the OP in a comment below).

NIntegrate[H[u]/u^3, {u, 0, 1}, MaxRecursion -> 12]
(* 21636.9 *)

For completeness, H (not H/u^3) is given by

Plot[H[u], {u, 0, 100}, PlotRange -> All, ImageSize -> Large, 
    AxesLabel -> {u, H}, LabelStyle -> Directive[12, Bold, Black]]

Second Addendum

A numerical solution can be obtained directly by substituting r[u] = u^2 w[u], as suggested by the symbolic solution above. The ODE becomes

(Unevaluated[H[u] == u*D[r[u], u] - 2 r[u]] /. r[u] -> u^2 w[u]) // Simplify
(* H[u] == u^3 Derivative[1][w][u] *)

As already noted above, r[10^6] == -1 is not a valid proxy for r[Infinity] == -1. So, arbitrarily choose w[1] == 1. Then, r[u] is given by

sol = u^2 w[u] /. Flatten@NDSolve[{w'[u] == H[u]/u^3, w[1] == 1}, w, {u, 10^-6, 2}, 
    AccuracyGoal -> Infinity]

plus C[1] u^2 as obtained symbolically earlier.

Plot[sol, {u, 10^-6, 2}, ImageSize -> Large, AxesLabel -> {u, r}, 
    LabelStyle -> Directive[12, Bold, Black]]