# To deal with infinity in NDSolve

Edit: According to Michael E2's comments, I did a few trail and error, and modify the question so that it uses exactly the same expressions I am having difficulties.

I just started to learn to use Mathematica, and was trying to obtain an inverse function. Helpful comments lead to the use of NDSolve, but I encountered the problem of infinity at the boundary of the integral.

f[r_] = (341/100 - (
92364544068700820988418540051678140032920439088266921797391767646\
86724351 (-1 + r)^14)/(
37439460288152388952575314975051235184793439401825645003319289114\
06568908800 r^14) - (
96562350197351134609086884969142247689186973725782763756859921174\
69 (-1 + r)^13)/(
60121690831554689932964817665258724812106472491722943248608203122\
49600 r^13) - (
133793861834595394886051762714719894875309418247714755031649667 \
(-1 + r)^12)/(
32499528808765052601032375726182804606730798541412729299084126720\
0 r^12) + (
8586490869886236032509596844465511159994386933570577541 (-1 +
r)^11)/(866273426372582901530614739357145019038988134825038599\
6800 r^11) + (
242824271326275467219778175017917876400335012677 (-1 + r)^10)/(
101170500983074227245353843599265287753122406259968 r^10) + (
38325063587426818959553212345333825733832573 (-1 + r)^9)/(
10875648319918368783953724555093403882956201600 r^9) + (
2707027650502829341747213307917371224783 (-1 + r)^8)/(
685879501776455635477799297139558154886400 r^8) + (
267199157832092112190906107808181 (-1 + r)^7)/(
89002913317143329270116576297209600 r^7) - (
14579284493690360915068655 (-1 + r)^6)/(
21048839588767223836466884944 r^6) - (
947101525623135764707 (-1 + r)^5)/(
84659318511493805075700 r^5) - (18202516748074031 (-1 + r)^4)/(
395488325322802800 r^4) - (9107664079253 (-1 + r)^3)/(
37412574526800 r^3) + (518274061 (-1 + r)^2)/(147465450 r^2) - (
4872 (-1 + r))/(775 r)) (-1 + r) r;

h[r_] = (341/100 - (
51413578086595983462378087137074096821605239174227636724824039221\
613454729 (-1 + r)^14)/(
41183406316967627847832846472556358703272783342008209503651218025\
472257996800 r^14) - (
13324089431333325317745097102721384152081575367154527677796135168\
9 (-1 + r)^13)/(
20613151142247322262730794628088705649865076282876437685237098213\
4272 r^13) + (
12238337278980590306571714477385268860836181893827453486014709 \
(-1 + r)^12)/(
64999057617530105202064751452365609213461597082825458598168253440\
r^12) + (
2054361172829064442171240720449805539268393215717818717 (-1 +
r)^11)/(173254685274516580306122947871429003807797626965007719\
9360 r^11) + (
3326409140357155826200623581130099710848389632453 (-1 + r)^10)/(
1517557514746113408680307653988979316296836093899520 r^10) + (
1583142347853489561731526395199754514115347 (-1 + r)^9)/(
543782415995918439197686227754670194147810080 r^9) + (
1900102987692718741210041373302547989269 (-1 + r)^8)/(
685879501776455635477799297139558154886400 r^8) + (
138947840710688646125560875614183 (-1 + r)^7)/(
267008739951429987810349728891628800 r^7) - (
3771160574785742691591658393 (-1 + r)^6)/(
526220989719180595911672123600 r^6) - (
5431784652331280941019 (-1 + r)^5)/(
169318637022987610151400 r^5) - (162489965633545399 (-1 + r)^4)/(
1186464975968408400 r^4) - (38932424921581 (-1 + r)^3)/(
37412574526800 r^3) + (411191516 (-1 + r)^2)/(73732725 r^2) - (
5847 (-1 + r))/(775 r)) (-1 + r) r;

MyPrecision = 500;
TortZero = 10^MyPrecision;
MyRecursion = 1000;


f[r_] is some complicated polynomial function, h[r_] is a similar one. Both functions satisfy f[1]=h[1]=1 and they approach $r^2$ at $r\rightarrow +\infty$. The function to be solved in question is $$r^*=r^*(r),$$ with $$\frac{dr^*}{dr}=\frac{1}{\sqrt{f(r)h(r)}},$$ satisfying the boundary condition $r^*(+\infty)=0$ so that $$r^*=(-1)\int_r^\infty dr'\frac{dr^*}{dr}(r').$$ It is noted that $\frac{dr^*}{dr}$ is singular at $r\rightarrow 1$ and $\frac{dr^*}{dr}\rightarrow \frac{1}{r^2}$ when $r\rightarrow +\infty$. Since the integral to infinity is convergent, it was possible to introduce the boundary condition as above. For a given $r$, the corresponding $r^*$ can be obtained by using numerical integral NIntegrate.

 rStar[r_?NumericQ] :=
NIntegrate[-1/(z^2 Sqrt[f[1/z] h[1/z]]), {z1, 0, 1/r}, WorkingPrecision -> MyPrecision,
MaxRecursion -> MyRecursion]


Though slow, it seems to work, Unfortunately, I need to obtain the inverse function at a very high precision. Following helpful suggestions/comments on this site, I tried to use NDSolve to deal with the problem with the following definitions

B2[z2_] = FullSimplify[-(z2)^2 Sqrt[f[1/z2] h[1/z2]], Assumptions -> {0 < z2 < 1/r0}] ;

C1[r1_] = FullSimplify[Sqrt[f[r1] h[r1]], Assumptions -> {0 < r1 < \[Infinity]}];


First I tried to solve the function r(r*) using the folllowing

s = NDSolve[{rI'[rS] == C1[rI[rS]], rI[-1/2] == 10},
rI[rS], {rS, -100, 0}, WorkingPrecision -> 100];


It did not work (even with the options used in the second case below), and then I thought maybe express the problem in terms of $z=\frac{1}{r}$ may help, and did a few trial and error. It seems to work with the following command.

s = NDSolve[{zI'[rS] == B2[zI[rS]], zI[-1/2] == 1/2},
zI[rS], {rS, -100, 0}, WorkingPrecision -> 500, AccuracyGoal -> 25,
PrecisionGoal -> 50, Method -> "Adams"];


My question is the following. I felt that most of the time I was trying without a well defined purpose, just trying out different options and different precision until I find a specific equation with specific settings that work. I also noticed that the version of Mathematica affect the results so that I eventually installed a old version, 9. Is there any guideline about what one should do when this sort of problem comes around? Many thanks in advance.

• Is there a mistake? The integral diverges (Integrate[1/Sqrt[f[r] h[r]], {r, 1, Infinity}]). Indeed, the integrand is asymptotic to a constant times 1/r (N@Limit[Sqrt[r^2/f[r]], r -> Infinity]). Jun 20, 2016 at 15:16
• Thanks for the answer. Yes, the integral is divergent if carried out from 1. But in practice, it is carried out from $r$ to $\infty$. I will edit the question to make this point clear. Jun 20, 2016 at 17:34
• I mean it's convergent at 1 and divergent at infinity. (Integrate[1/Sqrt[f[r] h[r]], {r, 10, Infinity}]) Jun 20, 2016 at 17:39
• I gave a fuller account of my thinking in an answer below. Jun 20, 2016 at 18:15
• I am really appreciated for the help. Clearly it was a mistake from my part. I edited the question, and will think about it according to your answer below. Jun 20, 2016 at 22:58

The singularity at r == 0 is of the type $O(1/\sqrt{r-1})$ and is easily removed with a substitution $r - 1 = u^2$. As for $r = \infty$, f[r] has a bunch of factors of the form $$(A + B r^n)/(C r^n)$$ plus the final two factors $(r-1)r$. This means that the function is asymptotic to $a r^2$ at $r = \infty$ for some constant $a$. This makes the integrand asymptotic to $b/r$, for $b = 1/\sqrt{a}$. Therefore the integral of $$\frac{dr^*}{dr}=\frac{1}{\sqrt{f(r)h(r)}}$$ will diverge (logarithmically) as $r \rightarrow \infty$, since we're taking $h \equiv 1$. Note that the general solution is $r^* = r^*_0(r) + C$, where $r^*_0$ is any particular solution; for example, it could be the solution to an IVP with a given value for $r^*_0(1)$. Then for the desired solution, we would pick $C$ so that $$\lim_{r \rightarrow \infty} r^*_0(r) + C = 0\,.$$ But this is impossible since $r^*_0(r)$ diverges at $r = \infty$.

One can see this computationally, too.

First, the substitution to get the new integrand:

int0 = 1/Sqrt@f[r] Dt[r] /. r -> u^2 + 1 /. Dt[u] -> 1 // Simplify[#, u > 0] &;


Then integrate:

{sol} = NDSolve[{rS1'[u] - int0 == 0, rS1[0] == 0}, rS1, {u, 0, 10^15},
"ExtrapolationHandler" -> {Nothing &, "WarningMessage" -> False},
WorkingPrecision -> 32, PrecisionGoal -> 10];


Various demonstrations of logarithmic growth:

Asymptotic growth rate: $\lim_{r \rightarrow \infty} {1 \over \sqrt{f(r)}}\big/{1 \over r} \approx 1.69859$.

coeff = N@Limit[Sqrt[r^2/f[r]], r -> Infinity]
(*  1.69859  *)


Comparison with the corresponding logarithm:

Plot[{rS1[Sqrt[r - 1]] /. sol, coeff*Log[r]}, {r, 1, 10^15}]


Plot[coeff*Log[r] - rS1[Sqrt[r - 1]] /. sol, {r, 1, 20},
PlotRange -> All, PlotLabel -> "Difference of rS1 and log"]


Constant rate of growth between successive powers of 10.

With[{valsAtPowersOfTen = rS1 /@ (10^Range[20]) /. sol // N},
Differences@valsAtPowersOfTen]
(*
{7.80958, 7.82216, 7.82229, 7.82229, 7.82229, 7.82229, 7.82229,
7.82229, 7.82229, 7.82229, 7.82229, 7.82229, 7.82229, 7.82229}
*)

• I am ashamed for my mistake. I thought to simplifies the problem but mistakenly ignores some important points. I will edit the question, and try to solve it to the best of my knowledge before asking for help. Jun 20, 2016 at 21:54
• @gamebm No problem. I've made such mistakes, too. Jun 21, 2016 at 0:12