Skip to main content
4 of 5
added second addendum
bbgodfrey
  • 62.1k
  • 18
  • 92
  • 160

The 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]]

enter image description here

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]]

enter image description here

Second Addendum

A numerical solution can be obtained directly by substituting r[u] = u^2 w[u], as suggested by the symbolic solution above. Then 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. The 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]]

enter image description here

bbgodfrey
  • 62.1k
  • 18
  • 92
  • 160