4
$\begingroup$

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.

$\endgroup$
6
  • $\begingroup$ 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]). $\endgroup$
    – Michael E2
    Jun 20, 2016 at 15:16
  • $\begingroup$ 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. $\endgroup$
    – gamebm
    Jun 20, 2016 at 17:34
  • $\begingroup$ I mean it's convergent at 1 and divergent at infinity. (Integrate[1/Sqrt[f[r] h[r]], {r, 10, Infinity}]) $\endgroup$
    – Michael E2
    Jun 20, 2016 at 17:39
  • $\begingroup$ I gave a fuller account of my thinking in an answer below. $\endgroup$
    – Michael E2
    Jun 20, 2016 at 18:15
  • $\begingroup$ 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. $\endgroup$
    – gamebm
    Jun 20, 2016 at 22:58

1 Answer 1

5
$\begingroup$

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

Mathematica graphics

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

Mathematica graphics

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}
*)
$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$
    – gamebm
    Jun 20, 2016 at 21:54
  • $\begingroup$ @gamebm No problem. I've made such mistakes, too. $\endgroup$
    – Michael E2
    Jun 21, 2016 at 0:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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