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


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

基于 $x_{g}$,评估函数 $(f^{'} \left ( x_{g} \right ))$ 的值。然后在该点绘制一条切线,使其与 𝑥 轴相交于 $x_{n}$。现在我们有两个点 $(x_{g},f\left ( x_{g} \right ))$ 和 $(x_{n} ,0)$。通过这些点经过的直线的斜率可以写成公式 1 所示。

$$\mathrm{f^{'} \left ( x_{g} \right )=\frac{0-f\left ( x_{g} \right )}{x_{n}-x_{g}}}$$

因此,$x_{n}$ 可以计算为 -

$$\mathrm{x_{n}-x_{_{g}}=-\frac{f\left ( x_{g} \right )}{f^{'}\left ( x_{g} \right )}}$$

$$\mathrm{x_{n}=x_{_{g}}-\frac{f\left ( x_{g} \right )}{f^{'}\left ( x_{g} \right )}}$$

现在,这个新的 x 值将作为下一步的猜测。基于这个新的猜测,再次评估函数的值,再次评估斜率,并重复该过程,直到获得根 $(i.e\left | x_{g} -x_{n}\right |<10^{-5})$。

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

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

假设我们想找到方程 $x^{2}+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)')

假设 𝑥 的初始猜测以开始第一次迭代。还将误差 $(\left | x_{g}-x_{n} \right |)$ 设置为大于收敛标准的值。在本文中,我将收敛标准设置为 $<10^{-5}$,但您可以根据所需的精度级别进行设置。并将循环计数器设置为 1。

# Initial Guess
xg=10

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

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

# 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

误差图如下所示 -

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

因此,对应于 $x_{g}=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

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

函数图将如下所示 -

因此,对应于 $x_{g}=-10$,根为 -3。

完整 Python 代码

完整代码如下 -

# 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

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

结论

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

更新时间: 2023年3月15日

2K+ 次浏览

开启您的 职业生涯

通过完成课程获得认证

开始学习
广告