In [1]:
import numpy as np
import matplotlib.pylab as plt

Solving the expression $1 - 2e^{-(R/h)} + R/h = 0.5$¶

First, let's simplify things with the substitution $x=R/h$ so that the equation becomes $1-2e^{-x}+x=0.5$

In [2]:
# note that 1 - 2e^(-x) + x = 0.5 is the same as 1 - 2e^(-x) + x - 0.5 = 0
# so we want to know the value of x where  1 - 2e^(-x) + x - 0.5 = 0.
# this is called "root-finding".

# define that function
def func(x):
    return 1 - 2*np.exp(-x) + x - 0.5

# and make a quick plot
x = np.arange(0,3,0.01)
plt.plot(x,func(x))
plt.gca().axhline(y=0,ls=':')
plt.xlabel('x')
plt.ylabel('f(x)')
Out[2]:
Text(0, 0.5, 'f(x)')

root-finding by determining the value in the array f(x) which is closest to 0¶

In [3]:
# evaluate t over a range of values that you think bracket the solution
# the step-size will determine the accuracy of your answer
x = np.arange(0,3,0.001) 

# find the index in the t array where func is closest to 0.
idx = np.argmin(np.abs(func(x))) 

# x[idx] is your solution
solution1=x[idx]
print('Solution is x = {:.3f}'.format(solution1)) 
Solution is x = 0.599

root-finding by using scipy's root-finding routines¶

In [4]:
from scipy.optimize import root

# you need to give the root-finder a rough guess at the correct value
sol_guess = 0.5  

# call the root-finding routine, printing out all its info
print('Printing out all the info from the root-finder:')
print('*****************************************')
print(root(func,sol_guess)) # give it the function and your initial guess for the solution
print('*****************************************')

# the solution is given in the x parameter of the root-finder, and we will
# take the first root value (which in this case is the only root value
solution2 = root(func,sol_guess).x[0]
print('Solution is x = {:.3f}'.format(solution2)) 
Printing out all the info from the root-finder:
*****************************************
    fjac: array([[-1.]])
     fun: array([-2.22044605e-16])
 message: 'The solution converged.'
    nfev: 7
     qtf: array([1.70319314e-12])
       r: array([-2.09886731])
  status: 1
 success: True
       x: array([0.59886728])
*****************************************
Solution is x = 0.599

So both ways give a solution of x=0.599.¶

Since $x=R/h$, this means that the solution is $R = 0.599h$¶