1
$\begingroup$

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!

$\endgroup$
  • $\begingroup$ The errors are from NIntegrate, not NDSolve. Also, I suggest you delete WorkingPrecision -> 12. $\endgroup$ – bbgodfrey Jul 27 '17 at 18:39
  • $\begingroup$ But when I just do up to plotting the NIntegrate part, I did get a graph with Only Warning that the integral is converging slowly. Not that it failed to converge. So I assume there must be something with the NDSolve. All the bad warnings pop up after NDSolve command... is there anything I can do? Thanks. $\endgroup$ – Duke Smith Jul 27 '17 at 19:54
  • $\begingroup$ Note: Above warning lines in the original post are for when I removed WorkingPrecision->12. When I put it in as the code That I wrote up there, it takes forever to run on my computer (cheap computer, possibly). I am a beginner in mathmatica $\endgroup$ – Duke Smith Jul 27 '17 at 20:04
1
$\begingroup$

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

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

enter image description here

$\endgroup$
  • $\begingroup$ Thank you! I've got a clue. It's kinda a different question, suppose I discard the u^2C[1] term, then i tried to NIntegrate H / u^3 term let's say 10^-9 to 100 (to avoid singularity at 0) it doesn't work too well. Do you know how to deal with this situation? I was expecting a convergent integral near 0 since H/u^3 decay slower than 1/u. $\endgroup$ – Duke Smith Jul 28 '17 at 2:29
  • $\begingroup$ @DukeSmith So that I am sure I understand, you would like to integrate H[u]/u^3 from zero to 100. Is that correct? $\endgroup$ – bbgodfrey Jul 28 '17 at 2:57
  • $\begingroup$ Yes, that is what I meant. u^2 * Integrate H[u]/u^3 Thank you for fast replies! $\endgroup$ – Duke Smith Jul 28 '17 at 3:01
  • $\begingroup$ @DukeSmith I do not have a quick answer for you now but will see what I can do tomorrow. $\endgroup$ – bbgodfrey Jul 28 '17 at 3:44
  • $\begingroup$ I followed your options for H[u], but I still get NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in u near {u} = {5.27879*10^-7}. NIntegrate obtained 1973.3583278644496` and 3.302641685197459` for the integral and error estimates. when I try to integrate for H[u]/u^3. Not sure how you did it. $\endgroup$ – Duke Smith Jul 29 '17 at 14:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.