Precision problem with numerical solution of a differential equation

Question

I want solve $$ 2\sqrt{|\gamma|}x = \int_{1}^{t} dy \sqrt{\frac{1+2|\gamma_2|y}{y^2(1-y)}} $$ where $0<t\leq $1. I'm using a for cycle to evaluate t, calculate the integral and the assign the $x$.

SfereDure = {}; 
For[t = 0.1, t <= 1, t = t + 0.1, 
   x = Integrate[Sqrt[(1 + 2*0.4*y)/(y^2*(1 - y))], {y, 1, t}]; 
 (*Print["t=",t]
 Print["x=", x]*)

 SfereDure = Prepend[SfereDure, {x/2, N[t, 1]}]; 
    SfereDure = Prepend[SfereDure, {-x/2, N[t, 1]}]]

I aspect $x=0$ at the last value of $t$ (i.e t=1), instead it isn't. I think that the problem is the precision of the calculus, in fact in the list "SfereDure"

{{-0.108449, 1.}, {0.108449, 1.}, {0.435952, 0.9}, {-0.435952, 0.9}, {0.635334, 0.8}, 
{-0.635334, 0.8}, {0.804738, 0.7}, {-0.804738, 0.7}, {0.965468, 0.6}, {-0.965468, 0.6}, 
{1.12849,0.5}, {-1.12849, 0.5}, {1.30375, 0.4}, {-1.30375, 0.4}, {1.50526, 0.3}, 
{-1.50526, 0.3}, {1.76102, 0.2}, {-1.76102, 0.2}, {2.15692, 0.1}, {-2.15692, 0.1}}

there is "1." that is 0.9999999999999999` and thus the integral $$ \int_{1}^{0.9999999999999999`}... $$ is not $0$.

Answer 1

I think your problem can be illustrated as simple as this:

res = {};
For[t = 0.1, t <= 1, t += 0.1, AppendTo[res, t]];
res // InputForm

and Mathematica is just doing what it is asked for: This is the most simple example of the ubiquitous problems with machine precision numbers (you might want to learn at least the very basics of doing numeric analysis with machine precision numbers). You might find it interesting that in Mathematica the following are simpler and in some sense do a better job (but of course can't avoid the machine precision problems completely):

Reap[Do[Sow[t], {t, 0.1, 1, 0.1}]] // InputForm

and the much appreciated:

Table[t, {t, 0.1, 1, 0.1}] // InputForm

If that still isn't good enough for your needs you have all the possibilities to get more precision:

Table[t, {t, 1/10, 1, 1/10}]
N[Table[t, {t, 1/10, 1, 1/10}]]
N[Table[t, {t, 1/10, 1, 1/10}],100] 

But be aware of the fact that higher precision (or exact solutions) usually come only at a price (usually computation time).

In case you are interested I have given some explanation here why I think that using For loops in Mathematica is almost never a good idea, your numerical problem being only one of the least important reasons...

You should also learn to use the numeric versions of the functions that Mathematica provides (here NIntegrate vs. Integrate). These are usually faster and work better when machine precision numbers are involved. The following seems to return what you expect, even though it works with just ordinary machine precision:

Reverse[Flatten[Table[
  x = NIntegrate[Sqrt[(1 + 2*0.4*y)/(y^2*(1 - y))], {y, 1, t}];
  {{x/2, t}, {-x/2, t}},
  {t, 0.1, 1, 0.1}
], 1]]

Answer 2

Your integral simply changes very rapidly near 1.

 Plot[ Evaluate[ 
         Integrate[Sqrt[(1 + 2 .4*y)/(y^2*(1 - y))], {y, t, 1}, 
         Assumptions -> {0 < t && t < 1}] ]  , {t, .5, 1}]

Your real issue I suppose is assuming your For loop ends at one. It does not due to numerical roundoff, the last value is ~ 1-10^-16 If you insist on using a For loop try:

 n = 10
 For[ it = 1 , it <= n , ++it ,
      t = N[it/n] ;
      x=Integrate...
   ]

Answer 3

Why, exactly, are you doing this? The integral can be evaluated in closed form:

Integrate[Sqrt[(1 + a y)/(y^2 (1 - y))], y]

Gives: $$ \frac{\sqrt{y-1} y \sqrt{\frac{a y+1}{y^2-y^3}} \left(\sqrt{a} \log \left(2 a y+2 \sqrt{a} \sqrt{y-1} \sqrt{a y+1}-a+1\right)-\tan ^{-1}\left(\frac{(a-1) y+2}{2 \sqrt{y-1} \sqrt{a y+1}}\right)\right)}{\sqrt{a y+1}} $$

And if you are opposed to symbolic methods, there is a perfectly good function NDSolve[]...

Answer 4

If you want better approximations, then do the integration with exact quantities and put off the numerical approximations as long as possible.

SfereDure = {};
For[t = 1/10, t <= 1, t += 1/10, 
  x = Integrate[Sqrt[(1 + 8/10*y)/(y^2*(1 - y))], {y, 1, t}];
  SfereDure = Prepend[SfereDure, {x/2, t} // N];
  SfereDure = Prepend[SfereDure, {-x/2, t} // N]]

SfereDure

{{0., 1.}, {0., 1.}, {0.435952, 0.9}, {-0.435952, 0.9}, 
 {0.635334, 0.8}, {-0.635334, 0.8}, {0.804738, 0.7}, {-0.804738, 0.7}, 
 {0.965468, 0.6}, {-0.965468, 0.6}, {1.12849, 0.5}, {-1.12849, 0.5}, 
 {1.30375, 0.4}, {-1.30375, 0.4}, {1.50526, 0.3}, {-1.50526, 0.3}, 
 {1.76102, 0.2}, {-1.76102, 0.2}, {2.15692, 0.1}, {-2.15692, 0.1}}