Error in newton raphson method finding root

Your main problem is to call something which is not callable.

• Functions are callable
• Sympy expressions are not callable

import sympy as sym

def f(x):
    return x**2

x = sym.Symbol("x")
g = 2*x

print(callable(f))  # True
print(callable(g))  # False
print(f(0))  # print(0)
print(g(0))  # error

So, in your case

def p(x):
    return 924*x**6 - 2772*x**5 + 3150*x**4 - 1680*x**3 + 420*x**2 - 42*x + 1

print(p(0))  # p is callable, gives the result 1
print(p(1))  # p is callable, gives the result 1
print(p(2))  # p is callable, gives the result 8989
print(callable(p))  # True

And now, if you use a symbolic variable from sympy you get:

x = sym.Symbol("x")
myp = p(x)
print(callable(myp))  # False
print(myp(0))  # gives an error, because myp is an expression, which is not callable

So, if you use diff to get the derivative of the function, you will get a sympy expression. You must transform this expression to a callable function. To do it, you can use the lambda function:

def p(x):
    return 924*x**6 - 2772*x**5 + 3150*x**4 - 1680*x**3 + 420*x**2 - 42*x + 1

x = sym.Symbol("x")
myp = p(x)
mydpdx = sym.diff(myp, x)  # get the derivative, once myp is an expression
dpdx = lambda x0: mydpdx.subs(x, x0)  # dpdx is callable
print(dpdx(0))  # gives '-42'

Another way to make a sympy expression callable is to use lambdify from sympy:

dpdx = sym.lambdify(x, mydpdx)

In sympy functions usually are represented as expressions. Therefore, sym.diff() returns an expression and not a callable function. It makes sense to also represent p as an expression involving x.

To "call" a function p on x, p.subs(x, value) is used. To get a numeric answer, use p.subs(x, value).evalf().

import sympy as sym

acc = 10 ** -4
x = sym.Symbol('x')
p = 924 * x ** 6 - 2772 * x ** 5 + 3150 * x ** 4 - 1680 * x ** 3 + 420 * x ** 2 - 42 * x + 1

p_prime = sym.diff(p)

def newton(p, p_prime, acc, start_val):
    xi = start_val
    delta = sym.oo
    while abs(delta) > acc:
        delta = (p / p_prime).subs(x, xi).evalf()
        print(xi, delta)
        xi -= delta
    return xi

a = newton(p, p_prime, 10**-4, 0)

Intermediate values:

0 -0.0238095238095238
0.0238095238095238 -0.00876459050881560
0.0325741143183394 -0.00117130331903862
0.0337454176373780 -1.98196463560775e-5
0.0337652372837341

PS: To plot a function represented as an expression:

sym.plot(p, (x, -0.03, 1.03), ylim=(-0.5, 1.5))

There are a few mistakes in your code, correcting them as pointed out in the following modified code snippet, it works fine:

def p_prime(xval):            # define as a function and substitute the value of x
    return sym.diff(p(x)).subs(x, xval)

def newton(p, p_prime, prec, start_val): # rename the output precision variable 
    x = start_val                        # otherwise it's shadowing acc tolerance
    delta = p(x) / p_prime(x)            # call the function p_prime()
    while abs(delta) > acc:
        delta = p(x) / p_prime(x)
        x = x - delta 
    return round(x, prec)                # return when the while loop terminates

a = newton(p, p_prime,4, 0.1)
a
# 0.0338

The following animation shows how Newton-Raphson converges to a root of the polynomial p(x).

Sign up or log in

Post as a guest