CODEWITHSUNDEEP
pythonnumpymathnewtons-method
snzm
snzm

Reputation: 159

Newton's method for finding roots

I'm trying to implement newton's method for finding roots in Python.

Expected result is point B, but instead Python returns point A:

geogebra

Code:

import matplotlib.pyplot as plt
import numpy as np

def f(theta):
    return 1 - (((2 * 1.5) * np.sin(theta))/ 2.7)

def derivative(f, x):
    dx = 1E-8
    return (f(x + dx) - f(x)) / dx

def x_next(f, x_n):
    return 1 - (f(x_n) / derivative(f, x_n))

def newtons_method(f, x_n = 1, i = 0, max_iter = 100):
    i = i + 1
    if (i == max_iter):
        return None
    x_n = x_next(f, x_n)
    if (abs(f(x_n)) < 1E-4):
        return x_n
    print("i:",i,"x_n:",x_n,"f(x_n)",f(x_n))
    newtons_method(f, x_n, i, max_iter)

print(newtons_method(f))

Upvotes: 0

Views: 1515

Answers (1)

Rory Daulton
Rory Daulton

Reputation: 22534

Your main problem is in your routine x_next. You have a 1 where there should be an x_n. So the routine should be

def x_next(f, x_n):
    return x_n - (f(x_n) / derivative(f, x_n))

Your derivative routine is also poor. If you have to approximate the derivative, Newton-Raphson is not the best method to use. Your approximation method that you use is also not good numerically, though it does follow the definition of derivative. If you must use an approximation, use

def derivative(f, x):
    dx = 1E-8
    return (f(x + dx) - f(x - dx)) / (2.0 * dx)

But in this case, the derivative is very easy to calculate directly. So it is better to use

def derivative(f, x):
    return -2 * 1.5 * np.cos(x) / 2.7

You also do not print your final approximation the the root and its function value--you calculate it and return without printing it. So place your print statement before you test it to return.

With those changes (plus commenting out the import of matplotlib which you never use), your code is now

#import matplotlib.pyplot as plt
import numpy as np

def f(theta):
    return 1 - (((2 * 1.5) * np.sin(theta))/ 2.7)

def derivative(f, x):
    return -2 * 1.5 * np.cos(x) / 2.7

def x_next(f, x_n):
    return x_n - (f(x_n) / derivative(f, x_n))

def newtons_method(f, x_n = 1, i = 0, max_iter = 100):
    i = i + 1
    if (i == max_iter):
        return None
    x_n = x_next(f, x_n)
    print("i:",i,"x_n:",x_n,"f(x_n)",f(x_n))
    if (abs(f(x_n)) < 1E-4):
        return x_n
    newtons_method(f, x_n, i, max_iter)

print(newtons_method(f))

and the result is only two lines

i: 1 x_n: 1.1083264212579311 f(x_n) 0.005607493777795347
i: 2 x_n: 1.1196379358595814 f(x_n) 6.373534192993802e-05

which is what you want. If you insist on using a numerical approximation for the derivative, use the version I give above and the result is slightly different:

i: 1 x_n: 1.10832642185337 f(x_n) 0.005607493482616466
i: 2 x_n: 1.1196379360265405 f(x_n) 6.373526104597182e-05

Upvotes: 2

Related Questions