Using Python To Interpret Results

I am trying to solve the following question. 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)-

Solution 1:

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()

Output

Using your code

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()

Vectorized approach

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()

Solution 2:

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

Your code is correct.

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

The limit of (sqrt(9+x)-3)/x for x→0 is 1/6 and for moderately small values of x you have different values that are approaching 1/6 BUT when x gets really small the numerator is affected by roundoff in the result of math.sqrt, hence the loss of convergence to the limit value.

For small values of x 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 x/6 is decreasing, ② there is a loss of precision in the calculation of the square root and ③ these small 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.

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,

Is there a better way to evaluate this expression?