牛顿拉夫森法求根时出错

Question

我尝试使用 newton raphson 方法计算函数的导数，但出现以下错误：

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)

错误是：

delta = p(x)/p_prime(x)
TypeError: 'Add' object is not callable

Answer 1

你的代码有一些错误，按照下面修改后的代码sn-p中指出的方法改正，效果很好：

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

以下动画展示了 Newton-Raphson 如何收敛到多项式 p(x) 的根。

Answer 2

在 sympy 中，函数通常表示为表达式。因此，sym.diff() 返回一个表达式而不是一个可调用函数。也可以将p 表示为涉及x 的表达式。

要在x 上“调用”函数p，使用p.subs(x, value)。要获得数字答案，请使用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)

中间值：

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

PS：绘制表示为表达式的函数：

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

Answer 3

您的主要问题是call 不是callable。

• 函数是可调用的
• Sympy 表达式不可调用

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

所以，在你的情况下

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

现在，如果您使用来自sympy 的symbolic variable，您会得到：

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

所以，如果你使用diff 得到函数的导数，你将得到一个sympy expression。您必须将此表达式转换为callable 函数。为此，您可以使用lambda 函数：

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'

使 sympy 表达式可调用的另一种方法是使用 sympy 中的lambdify：

dpdx = sym.lambdify(x, mydpdx)