使用 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 中绘制切线并显示根的收敛。