# Precision problem with numerical solution of a differential equation

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

• What are $|\gamma|$ and $|\gamma_2|$? Are you going to solve it with respect to ...??? One equation, two unknowns? – Artes Sep 12 '14 at 14:27
• The mean of $|\gamma|$ and $\gamma_2$ is no important... The problem is the precision of the calculus – apt45 Sep 12 '14 at 14:31
• It goes to 0 for me if I increase the number of 9 in 0.99999. Which version of Mathematica do you use? May be you want to use NItegrate if you don't need an exact formula? – ybeltukov Sep 12 '14 at 14:54
• If you see the code, the variable increase of 0.1 at each step. So i expect that the last value of $t$ is 1, instead of $0.99999999999'$. I use Mathematica 9. What do you mean "if i increase the number of 9 in 0.9999"? – apt45 Sep 12 '14 at 14:57

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


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


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

• i'm doing this because i'm not interested in a exact solution. – apt45 Sep 12 '14 at 14:30
• What ARE you interested in? – Igor Rivin Sep 12 '14 at 14:34
• I don't understand why the for cycle doesn't work well, why the integral isn't 0 when $t=1$. – apt45 Sep 12 '14 at 14:36

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