Homework 9

We are going to explore the zeros of the zeroth bessel function of first kind, by using the bisection and Newton-Raphson methods and then plotting the comparison between them.

First we import the libraries and save the function on a certain interval just for a few points,

In [9]:
import matplotlib.pylab as plt
import numpy as np
from scipy import special
x=np.linspace(0,10,10)
y=special.j0(x)

So the plot goes as

In [10]:
fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(x,y,'.',label='data')
ax.set_title("Bessel Function $j_0(x)$")
ax.set_xlabel("$x$")
plt.legend()
plt.show()

To find the zeros, we define the function and its first derivative

In [11]:
def func(x):
    return special.j0(x)
def dfunc(x):
    return -special.j1(x)

Exploring the data, we define 3 intervals where the zeros are,

In [12]:
intervals=[[1,3],[4.5,6.5],[8,10]]

Then, we define the bisection and Newton-Raphson methods to return lists of the 10 firsts values,

In [13]:
def bisec(f,an,bn):
    pn=[]
    for i in range(10):
        pm=(an+bn)/2.
        if (f(an)*f(pm)>0):
            an=pm
        else:
            bn=pm
        pn.append(pm)
    return pn
def Newton(f0,df0,xi):
    xn=[]
    for i in range(10):
        xi-=f0(xi)/df0(xi)
        xn.append(xi)
    return xn

We Use both methods for the three intervals and plot the solutions.

In [14]:
final_newton=[]
final_bisection=[]
for i in intervals:
    xmin=i[0]
    xmax=i[1]
    data=bisec(func,xmin,xmax)
    data2=Newton(func,dfunc,xmin)
    plt.plot(data,'.-',label='Bisection')
    plt.plot(data2,'.-',label='Newton')
    plt.xlabel('Iteration')
    plt.ylabel('Result')
    plt.legend()
    plt.show()
    final_newton.append(data2[-1])
    final_bisection.append(data[-1])

Conclusions

Plotting this, we can clearly see that, both methods converge to the same result, but the Newton's method does it way faster.

To test our solutions, we can use the SciPy solutions as if they were the real data.

In [7]:
real=special.jn_zeros(0,3)
bis=np.array(final_bisection)
new=np.array(final_newton)

So, one can see the percentual error, so that

In [8]:
print((new-real)/real)
print((bis-real)/real)

[0. 0. 0.]
[-2.19842431e-04  2.54754495e-04  6.57476286e-05]

And as you can see, the Newton's method has exact the same result!.