Question:

I am trying to solve the following question.

<a href="https://i.stack.imgur.com/J0q5Y.jpg" rel="nofollow"><img alt="enter image description here" class="b-lazy" data-src="https://i.stack.imgur.com/J0q5Y.jpg" data-original="https://i.stack.imgur.com/J0q5Y.jpg" src="https://etrip.eimg.top/images/2019/05/07/timg.gif" /></a>

Here is my code I have produced

```
import numpy as np
import math
sum = 4
while sum <= 13:
b = 10**(-sum)
x = (math.sqrt(9+b)-3)/(b)
print(x)
sum +=1
```

which gives the following results

```
0.16666620370475727
0.1666666203714584
0.1666666618049817
0.1666666671340522
0.1666666804567285
0.1666666804567285
0.1666666804567285
0.1666666804567285
0.16653345369377348
0.16431300764452317
```

I am not sure if my code is correct or not. When I plug in 13 for `n`

into the original equation in wolfram I am getting something different. I figured as it got closer to 13 it would approach 0.1666666.

Also how would I make a graph from this? I guess this would be the best way to observe my results.

Answer1:Here is a complete solution together with the plot. Explanation: `np.logspace(4, 13, 10)`

creates the values of `x`

as 10^(4), 10^(5), 10^(6)... 10^(13). You take the inverse of its output to get your desired x-points as 10^(-4), 10^(-5), 10^(-6)... 10^(-13). You then loop over these x-values and solve for your equation. You save the output value for each x in a list and then plot it.

There are other vectorized approaches without having to create a loop. But this should get you started.

```
import numpy as np
import math
import matplotlib.pyplot as plt
xmesh = 1./np.logspace(4, 13, 10)
result = []
for x in xmesh:
elem = (math.sqrt(9+x)-3)/(x) # removed 'b' because xmesh is already inverse
result.append(elem) # Adds each element to the list for later plotting
plt.semilogx(xmesh, result, '-bo') # Uses a logarithmic x-axis
plt.gca().invert_xaxis() # Inverts the x-axis because you want it so as per comments
plt.xlabel('1/x', fontsize=16)
plt.ylabel('y', fontsize=16)
plt.show()
```

<strong>Output</strong>

<a href="https://i.stack.imgur.com/ePdq2.png" rel="nofollow"><img alt="enter image description here" class="b-lazy" data-src="https://i.stack.imgur.com/ePdq2.png" data-original="https://i.stack.imgur.com/ePdq2.png" src="https://etrip.eimg.top/images/2019/05/07/timg.gif" /></a>

<strong>Using your code</strong>

Following is how to make it work using your code without much modification

```
import numpy as np
import math
import matplotlib.pyplot as plt
sum = 4
xmesh = 1./np.logspace(4, 13, 10)
result = []
while sum <= 13:
b = 10**(-sum)
x = (math.sqrt(9+b)-3)/(b)
result.append(x)
sum +=1
plt.semilogx(xmesh, result, '-bo')
plt.gca().invert_xaxis()
plt.xlabel('1/x', fontsize=16)
plt.ylabel('y', fontsize=16)
plt.show()
```

<strong>Vectorized approach</strong>

```
import numpy as np
xmesh = 1./np.logspace(4, 13, 10)
result = (np.sqrt(9+xmesh)-3)/(xmesh) # No need to loop. Used `np.sqrt`
plt.semilogx(xmesh, result, '-bo')
plt.gca().invert_xaxis()
```

Answer2:## <blockquote>

I am not sure if my code is correct or not.

</blockquote><strong>Your code is correct.</strong>

What you see is the result of finite precision of floating numbers.

The limit of <em>`(sqrt(9+x)-3)/x`

</em> for <em>`x→0`

</em> is <em>`1/6`

</em> and for moderately small values of <em>`x`

</em> you have different values that are approaching <em>`1/6`

</em> BUT
when <em>`x`

</em> gets really small the numerator is affected by <em>roundoff</em> in the result of `math.sqrt`

, hence the loss of convergence to the limit value.

For small values of <em>x</em> the numerator can be approximated using the binomial approximation, (9+x)^0.5-3= 3*(1+x/9)^0.5-3 ≈ 3*(1+x/18)-3 = 3+x/6-3 = x/6, let's see, using Python, what is really happening

```
>>> for n in range(7,16):
... x = 10**(-n)
... print('%2d %10e %s %15e %15e'%(n, x, repr(sqrt(9+x)), sqrt(9+x)-3, x/6))
...
7 1.000000e-07 3.0000000166666667 1.666667e-08 1.666667e-08
8 1.000000e-08 3.000000001666667 1.666667e-09 1.666667e-09
9 1.000000e-09 3.0000000001666667 1.666667e-10 1.666667e-10
10 1.000000e-10 3.0000000000166667 1.666667e-11 1.666667e-11
11 1.000000e-11 3.0000000000016667 1.666667e-12 1.666667e-12
12 1.000000e-12 3.0000000000001665 1.665335e-13 1.666667e-13
13 1.000000e-13 3.0000000000000164 1.643130e-14 1.666667e-14
14 1.000000e-14 3.0000000000000018 1.776357e-15 1.666667e-15
15 1.000000e-15 3.0000000000000004 4.440892e-16 1.666667e-16
```

I've used `repr(...)`

to show the significant digits in the square root result and you can see different concurrent phenomena ① the number of significant digits used to represent <em>x</em>/6 is decreasing, ② there is a loss of precision in the calculation of the square root and ③ these <em>small</em> effects are amplified by a cancellation effect, as you subtract 3 from a substantially correct result `(3+δ)`

.

Also how would I make a graph from this? I guess this would be the best way to observe my results.

</blockquote>As you can see, it is not necessary to graph the results to gain an insight on the problem!!!

What remains to do, when you have recognized the roundoff caused instability, is to tackle the last part of your question,

<blockquote> <blockquote>Is there a better way to evaluate this expression?

</blockquote> </blockquote>