使用 Python 建模牛顿-拉夫森方法

在本教程中，我将向您展示如何借助称为牛顿-拉夫森方法的数值方法来评估多项式或超越方程的根。这是一种迭代方法，我们从一个初始猜测（自变量）开始，然后根据猜测评估 𝑥 的新值。并且该过程持续进行，直到达到收敛。该方法通过下图所示的图表进行解释。

基于 xg，评估函数 (f′(xg)) 的值。然后在该点绘制一条切线，使其与 𝑥 轴相交于 xn。现在我们有两个点 (xg,f(xg)) 和 (xn,0)。通过这些点经过的直线的斜率可以写成公式 1 所示。

f′(xg)=0−f(xg)xn−xg

因此，xn 可以计算为 -

xn−xg=−f(xg)f′(xg)

xn=xg−f(xg)f′(xg)

现在，这个新的 x 值将作为下一步的猜测。基于这个新的猜测，再次评估函数的值，再次评估斜率，并重复该过程，直到获得根 (i.e|xg−xn|<10−5)。

该方法速度很快，但一次只能给你一个根。要获得另一个根，您必须从另一个猜测开始并再次重复该过程。

牛顿-拉夫森方法的 Python 实现

假设我们想找到方程 x2+4x+3=0 的根。牛顿-拉夫森方法的 Python 实现如下 -

导入包 -

from pylab import *

只使用了一个模块，即 pylab，因为它包含 numpy。因此，无需单独导入它。

形成多项式及其导数函数，即 𝑓(𝑥) 和 𝑓'(x)。

f=lambda x: x**2+4*x+3
dfdx=lambda x: 2*x+4

我使用了 'lambda'，因为函数中只有一条语句。如果需要，也可以使用 'def' 方法。

使用 "linspace" 函数为 "x" 创建一个数组。

# Array of x
x=linspace(-15,10,50)

现在，此步骤是可选的。考虑适当的域绘制函数。我还将向您展示如何绘制切线，以及解决方案如何收敛。因此，如果您对视觉效果感兴趣，则可以遵循此步骤。

# Plotting the function
figure(1,dpi=150)
plot(x,f(x))
plot([-15,10],[0,0],'k--')
xlabel('x')
ylabel('f(x)')

假设 𝑥 的初始猜测以开始第一次迭代。还将误差 (|xg−xn|) 设置为大于收敛标准的值。在本文中，我将收敛标准设置为 <10−5，但您可以根据所需的精度级别进行设置。并将循环计数器设置为 1。

# Initial Guess
xg=10

# Setting initial error and loop counter
error=1
count=1

在 "for" 循环中，使用上述收敛标准求解公式 (2)。此外，绘制误差和切线。切线绘制在名为 figure(1) 的绘图中，误差在 figure(2) 中。此外，还将 xg 和 f(xg) 制成表格以显示不同的值。

# For printing x_g and f(x_g) at different steps
print(f"{'xg':^15}{'f(xg)':^15}")
print("===========================")

# Starting iterations
while error>1.E-5:
   # Solving Eq. 1
   xn=xg-f(xg)/dfdx(xg)
    
   # Printing x_g and f(x_g)
   print(f'{round(xg,5):^15}{round(f(xg),5):^15}')
   
   # Plotting tangents
   figure(1,dpi=300)
   plot([xg,xn],[f(xg),0])
   plot([xn,xn],[0,f(xn)],'--',label=f'xg={xg}')
   legend(bbox_to_anchor=(1.01, 1),loc='upper left', borderaxespad=0)
   
   # Evaluating error and plotting it
   error=abs(xn-xg)
   figure(2,dpi=300)
   semilogy(count,error,'ko')
   xlabel('Iteration count')
   ylabel('Error')
   
   # Setting up new value as guess for next step
   xg=xn
   
   # Incrementing the loop counter
   count=count+1
# printing the final answer
print("===========================")
print("\nRoot ",round(xn,5))
show()

收敛后，打印根。并显示绘图。

在上述情况下，我将初始猜测设为 10。因此，程序输出将如下所示 -

      xg                    f(xg)
======================================
      10                     143
    4.04167             35.50174
    1.10359              8.63228
    -0.2871               1.93403
   -0.85165              0.31871
   -0.99042              0.01926
   -0.99995               9e-05
     -1.0                    0.0
========================================
Root -1.0

误差图如下所示 -

带有切线的函数图在下图中显示。

因此，对应于 xg=10，根为 -1。对于第二个根，我们必须更改猜测，将其设为 -10。然后程序输出将如下所示 -

      xg          f(xg)
===========================
       -10          63
     -6.0625     15.50391
    -4.15433      3.64112
    -3.30925      0.71415
    -3.03652      0.07438
    -3.00064      0.00129
    -3.0            0.0
===========================
Root -3.0

现在，误差图将如下所示 -

函数图将如下所示 -

因此，对应于 xg=−10，根为 -3。

完整 Python 代码

完整代码如下 -

Open Compiler
# Importing module
from pylab import *

# Funciton for f(x) and f'(x)
f = lambda x: x ** 2 + 4 * x + 3
dfdx = lambda x: 2 * x + 4

# Array of x
x = linspace(-15, 10, 50)

# Plotting the function
figure(1, figsize=(7.20, 4.00))
plot(x, f(x))
plot([-15, 10], [0, 0], 'k--')
xlabel('x')
ylabel('f(x)')

# Initial Guess
xg = 10

# Setting initial error and loop counter
error = 1
count = 1

# For printing x_g and f(x_g) at different steps
print(f"{'xg':^15}{'f(xg)':^15}")
print("===========================")

# Starting iterations
while error > 1.E-5:
   # Solving Eq. 1
   xn = xg - f(xg) / dfdx(xg)

   # Printing x_g and f(x_g)
   print(f'{round(xg, 5):^15}{round(f(xg), 5):^15}')

   # Plotting tangents
   figure(1)
   plot([xg, xn], [f(xg), 0])
   plot([xn, xn], [0, f(xn)], '--', label=f'xg={xg}')
   legend(bbox_to_anchor=(0.4, 1.1), loc='upper left', borderaxespad=0)

   # Evaluating error and plotting it
   error = abs(xn - xg)
   figure(2, figsize=(7.20, 4.00))
   semilogy(count, error, 'ko')
   xlabel('Iteration count')
   ylabel('Error')

   # Settingup new value as guess for next step
   xg = xn

   # Incrementing the loop counter
   count = count + 1

# printing the final answer
print("===========================")
print("\nRoot ", round(xn, 5))
show()

您可以将代码直接复制到 Jupyter Notebook 中并运行它。

对于您选择的多项式，您可以更改上述代码中所示的函数和导数多项式，并根据您的猜测值，您将获得输出。例如，如果您想找到 𝑥3−sin2(𝑥)−𝑥=0 的根，则在上述代码中，函数及其导数将更改为 -

# Function for f(x) and f'(x)
f=lambda x: x**3-(sin(x))**2-x
dfdx=lambda x: 3*x**2-2*sin(x)*cos(x)-1

然后对于 1 的猜测值，程序输出将为 -

       xg           f(xg)
===========================
        1          -0.70807
     1.64919        1.84246
     1.39734        0.36083
     1.31747        0.0321
     1.30884        0.00036
     1.30874          0.0
===========================
Root 1.30874

并且，函数图将如下所示 -

Learn Python in-depth with real-world projects through our Python certification course. Enroll and become a certified expert to boost your career.

结论

在本教程中，您学习了如何使用牛顿-拉夫森方法求解方程的根。您还学习了如何在 pyplots 中绘制切线并显示根的收敛。

更新时间： 2023年3月15日

2K+ 次浏览

开启您的 职业生涯

通过完成课程获得认证

开始学习
广告