1

I was trying to use the newton raphson method to compute the derivative of a function and I got the following error:

import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
acc = 10**-4

x = sym.Symbol('x')
def p(x): #define the function
    return 924*x**6 - 2772*x**5 + 3150*x**4 - 1680*x**3 + 420*x**2 - 42*x + 1

p_prime = sym.diff(p(x))
def newton(p,p_prime, acc, start_val):
    x= start_val
    delta = p(x)/p_prime(x)
    while abs(delta) > acc:
        delta = p(x)/p_prime(x)
        print(abs(delta))
        x =x- delta;
        return round(x, acc);

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

The error was:

delta = p(x)/p_prime(x)
TypeError: 'Add' object is not callable
Sandipan Dey
  • 21,482
  • 2
  • 51
  • 63

3 Answers3

1

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

enter image description here

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

Sandipan Dey
  • 21,482
  • 2
  • 51
  • 63
0

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)
Carlos Adir
  • 452
  • 3
  • 9
0

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

plot of symbolic function

JohanC
  • 71,591
  • 8
  • 33
  • 66
  • Thanks, it worked. Can you elaborate what delta = sym.oo does? – The Emerging Star Dec 04 '21 at 13:55
  • And also how to display the grids, legends in the plot? – The Emerging Star Dec 04 '21 at 14:35
  • `oo` is infinity in sympy. Initializing delta with a very high value avoids having to put the calculation of delta twice. "Don't Repeat Yourself" (DRY) is an important principle. If the way to calculate ia changed in the future, and the calculation appears twice, it is very common to only change it at one spot or to make another mistake. – JohanC Dec 04 '21 at 15:35
  • Sympy's plotting is quite limited. It would be easier to use lambdify and plot with matplotlib and numpy. – JohanC Dec 04 '21 at 15:39