通过Python实现5个常用数值近似算法

二分法

二分法是一种简单有效的数值型迭代算法，对于一个在区间$\left[a,b\right]$上的连续函数fx，若满足$f(a)\cdot f(b)<0$，那么fx在$\left[a,b\right]$上必有根，此时设定分割点$x_0=\left(a+b\right)/2$将区间等分为$\left[a,x_0\right]$和$\left[x_0,b\right]$，计算$f\left(a\right)\cdot f\left(x_0\right)$、$f\left(x_0\right)\cdot f\left(b\right)$，则根必位于其中一区间内，重复上述过程N次直至根所在的区间长度足够小时停止迭代。尽管二分法能够保证最终收敛到根附近，但它也存在着一个极端情况，即每次迭代中根都恰好位于另一区间内，因此时间复杂度为$O={\log }_2\left(N\right)$。二分法的流程如下图所示

• 源码

import math
import numpy as np
import matplotlib.pyplot as plt

x = np.arange(-0.5, 4, 0.01)
y = np.cos(x)
plt.figure(figsize=(10, 8))
plt.plot(x, y,label="f(x)",linewidth=2)
plt.vlines(x=0, ymin=-1.5, ymax=1.5)
plt.hlines(y=0, xmin=-1, xmax=5)
#plt.xlabel("x")
#plt.ylabel("cos(x)")
plt.vlines(x=0.3, ymin=0, ymax=cos(0.3), linestyle="--", color="black")
plt.scatter(0.3, cos(0.3), c = 'yellow', marker = 'o',s=100, edgecolors="black",alpha = 1, linewidths = 1.5)
plt.vlines(x=3.5, ymin=0, ymax=cos(3.5), linestyle="--", color="black")
plt.scatter(3.5, cos(3.5), c = 'yellow', marker = 'o',s=100, edgecolors="black",alpha = 1, linewidths = 1.5)
plt.scatter(0.27, -0.05, c = 'black', marker = '$a$', s=80, alpha = 1, linewidths = 0.5)
plt.scatter(3.5, 0.08, c = 'black', marker = '$b$', s=120, alpha = 1, linewidths = 0.5)
plt.scatter(0.27, cos(0.3)+0.1, c = 'black', marker = '$f(a)$', s=580, alpha = 1, linewidths = 0.5)
plt.scatter(3.5, cos(3.5)-0.1, c = 'black', marker = '$f(b)$', s=580, alpha = 1, linewidths = 0.5)
plt.scatter(np.pi/2-0.05, cos(np.pi/2)-0.1, c = 'black', marker = '$x_0$', s=280, alpha = 1, linewidths = 0.5)
plt.scatter(np.pi/2, cos(np.pi/2), c = 'red', marker = 'o', s=200, edgecolors="black",alpha = 1,linewidths = 1.5)

plt.scatter((0.3+3.5)/2, 0, c = 'red', marker = 'o', s=100, edgecolors="black", alpha = 1, linewidths = 1.5)
plt.scatter((0.3+3.5)/2, 0+0.1, c = 'red', marker = '$a_1$', s=150, edgecolors="black", alpha = 1, linewidths = 0.5)
plt.scatter((0.3+3.5)/2-0.1, cos((0.3+3.5)/2)-0.1, c = 'black', marker = '$f(a_1)$', s=600, edgecolors="black", alpha = 1, linewidths = 0.5)
plt.scatter((0.3+3.5)/2, cos((0.3+3.5)/2), c = 'yellow', marker = 'o', s=100, edgecolors="black", alpha = 1, linewidths = 1.5)
plt.vlines(x=(0.3+3.5)/2, ymin=0, ymax=cos((0.3+3.5)/2), linestyle="--",color="black")
plt.scatter((0.3+3.5)/4, 0, c = 'red', marker = 'o', s=100, edgecolors="black", alpha = 1, linewidths = 1.5)
plt.scatter((0.3+3.5)/4, 0-0.1, c = 'red', marker = '$a_2$', s=150, edgecolors="black", alpha = 1, linewidths = 0.5)
plt.vlines(x=(0.3+3.5)/4, ymin=0, ymax=cos((0.3+3.5)/4), linestyle="--", color="black")
plt.scatter((0.3+3.5)/4, cos((0.3+3.5)/4), c = 'yellow', marker = 'o', s=100, edgecolors="black", alpha = 1, linewidths = 1.5)
plt.scatter((0.3+3.5)/4+0.05, cos((0.3+3.5)/4)+0.1, c = 'black', marker = '$f(a_2)$', s=600, edgecolors="black", alpha = 1, linewidths = 0.5)
plt.scatter(1.5, 0, c='red', marker = 'o', s=100, edgecolors="black", alpha = 1, linewidths = 1.5)
plt.scatter(1.38, 0.06, c='black', marker = '$a_3$', s=150, edgecolors="black", alpha = 1, linewidths = 0.5)
plt.scatter(-0.1, -0.1, c='black', marker = '$0$', s=150, edgecolors="black", alpha = 1, linewidths = 0.5)
plt.scatter(5.05, 0, c='black', marker = '>', s=100, edgecolors="black", alpha = 1, linewidths = 0.5)
plt.scatter(0, 1.53, c='black', marker = '^', s=100, edgecolors="black", alpha = 1, linewidths = 0.5)
plt.scatter(-0.2, 1.62, c='black', marker = '$f(x)$', s=600, edgecolors="black", alpha = 1, linewidths = 0.5)
plt.scatter(5.0, -0.09, c='black', marker = '$x$', s=130, edgecolors="black", alpha = 1, linewidths = 0.5)
plt.yticks(fontstyle = "normal", fontproperties = 'Times New Roman', size = 13)
plt.xticks(fontstyle = "normal", fontproperties = 'Times New Roman', size = 13)
#plt.xlabel('x', fontdict=fontdict_ita)
#plt.ylabel("f (x)", fontdict=fontdict_ita)
plt.legend(prop=fontdict_ita)
plt.tight_layout()
plt.show()

牛顿迭代法

牛顿拉夫逊迭代算法又简称为牛顿法，是我们最为熟悉且应用范围最广的数值型计算方法。牛顿迭代法是求方程根的重要方法之一，其最大优点是在方程$f\left(x\right)=0$的单根附近具有平方收敛，而且该算法还可以用来求方程的重根、复根等，此时线性收敛，但是可通过一些方法变成超线性收敛。另外该方法广泛用于计算机编程中。若设鞍点方程为$f\left(t\right)$，$r$是该方程的根，选择$x_0$作为迭代初始值，过点$\left(x_0,f\left(x_0\right)\right)$做曲线$f\left(t\right)$的切线$L$，则根据点斜式有L的表达式为
$$L:y=f(x_0)+f^{\prime }(x_0)(x-x_0)$$
该切线与x轴交于
$$x_1=x_0-\frac{f(x_0)}{f^{\prime }(x_0)}$$
过点$\left(x_1,f\left(x_1\right)\right)$再做曲线$f\left(x\right)$的切线得$x_2$，重复上述过程$N$次直至$x_{n+1}-x_n\le \epsilon$结束迭代。已经有大量的文献可以证明若方程的根孤立存在，则根的附近存在一个邻域，位于该邻域内的初始迭代值最终必收敛至方程的根。下图给出了牛顿法计算过程的示意

x = np.arange(-1, 5, 0.01)
k = 0.55
def f(x):
    return np.exp(k*x)-2.5
def d_f(x):
    return k*exp(k*x)
x0 = 4.5

plt.figure(figsize=(10, 8))
plt.plot(x, f(x), label="f(x)", linewidth=2)
plt.vlines(x=0, ymin=-5, ymax=15, color="black")
plt.hlines(y=0, xmin=-1, xmax=5, color="black")
n=5
for i in range(1, n):
    #cut = -f(x1)*(x1 - x0)/(f(x1) - f(x0)) + x1
    x1 = x0 - f(x0)/d_f(x0)
    xx = np.arange(x1, 5, 0.01)
    yy = f(x0) + d_f(x0)*(xx - x0)
    plt.plot(xx, yy, color="black", linestyle="--", linewidth=0.8)
    plt.scatter(x0, f(x0), c='yellow', marker = 'o', s=size_o , edgecolors="black", alpha = 1, linewidths = 1.5)
    plt.scatter(x0 - 0.18, f(x0) + 0.8, c='black', marker=marker_fx[i-1], s=800 , edgecolors="black", alpha = 1, linewidths = 0.5)
    plt.scatter(x0, -0.5, c='black', marker = marker_x[i-1], s=size_x , edgecolors="black", alpha = 1, linewidths = 0.5)
    plt.vlines(x=x0, ymin=0, ymax=f(x0), linestyle="--")
    x0 = x1
    x1 = cut

plt.yticks(fontstyle = "normal", fontproperties = 'Times New Roman', size = 13)
plt.xticks(fontstyle = "normal", fontproperties = 'Times New Roman', size = 13)
plt.scatter(5.05, 0, c='black', marker = '>', s=100, edgecolors="black", alpha = 1, linewidths = 0.5)
plt.scatter(0, 15.3, c='black', marker = '^', s=100, edgecolors="black", alpha = 1, linewidths = 0.5)
plt.scatter(-0.2, 16.2, c='black', marker = '$f(x)$', s=600, edgecolors="black", alpha = 1, linewidths = 0.5)
plt.scatter(5.0, -0.7, c='black', marker = '$x$', s=130, edgecolors="black", alpha = 1, linewidths = 0.5)
plt.legend(prop=fontdict_ita)
plt.tight_layout()
plt.show()

割线法

割线法是牛顿法的一种改进，其思想是将牛顿法中所求方程的切线替换为割线，再将割线与x轴的交点作为根的近似解。给定函数$f\left(x\right)$和函数上两点$\left(x_1,f\left(x_1\right)\right)$、$\left(x_2,f\left(x_2\right)\right)$，则割线L由两点式给出
$$y-f(x_1)=\frac{f(x_1)-f(x_0)}{x_1-x_0}(x-x_1)$$
此时L与$x$轴交于$x_2$点，重新做过$x_1$、$x_2$的割线与横轴交于$x_3$，重复上述过程$N$次直至$x_n-x_{n-1}$收敛于给定的误差限。与牛顿法不同的是，割线法需要预先给定两个初始迭代点。其算法过程如下图所示

x = np.arange(-1, 5, 0.01)
k = 0.7

def f(x):
    return np.exp(k*x)-2.5
def d_f(x):
    return k * exp(k * x)


x0 = 4.8; f0 = f(x0)
x1 = 3.8; f1 = f(x1)
n = 7

plt.figure(figsize=(10, 8))
plt.plot(x, f(x), linewidth=2.5, label="f(x)")
plt.vlines(x=0, ymin=-5, ymax=f(max(x))+1, color="black")
plt.hlines(y=0, xmin=-1, xmax=5, color="black")

marker_fx = []
marker_x = []
for i in range(1, n):
    lab = "$f(x_" + str(i) + ")$"
    lab2 = "$x_" + str(i) + "$"
    marker_fx.append(lab)
    marker_x.append(lab2)

for i in range(1, n):
    cut = -f(x1)*(x1 - x0)/(f(x1) - f(x0)) + x1
    xx = np.arange(cut, 5, 0.01)
    yy = (xx - x1)*((f(x1) - f(x0))/(x1 - x0)) + f(x1)
    plt.plot(xx, yy, color="black", linestyle="--", linewidth=0.8)
    plt.scatter(x0, f(x0), c='yellow', marker = 'o', s=size_o , edgecolors="black", alpha = 1, linewidths = 1.5)
    plt.scatter(x0 - 0.25, f(x0) + 1, c='black', marker=marker_fx[i-1], s=700 , edgecolors="black", alpha = 1, linewidths = 1)
    plt.scatter(x0, -1.0, c='black', marker = marker_x[i-1], s=size_x , edgecolors="black", alpha = 1, linewidths = 0.5)
    plt.vlines(x=x0, ymin=0, ymax=f(x0), linestyle="--")
    x0 = x1
    x1 = cut
    
plt.yticks(fontstyle = "normal", fontproperties = 'Times New Roman', size = 13)
plt.xticks(fontstyle = "normal", fontproperties = 'Times New Roman', size = 13)
plt.scatter(5.05, 0, c='black', marker = '>', s=100, edgecolors="black", alpha = 1, linewidths = 0.5)
plt.scatter(0, 31.3, c='black', marker = '^', s=100, edgecolors="black", alpha = 1, linewidths = 0.5)
plt.scatter(-0.2, 33.2, c='black', marker = '$f(x)$', s=600, edgecolors="black", alpha = 1, linewidths = 0.5)
plt.scatter(5.0, -1.2, c='black', marker = '$x$', s=130, edgecolors="black", alpha = 1, linewidths = 0.5)
plt.legend(prop=fontdict_ita)
plt.tight_layout()
plt.show()

弦截法

弦截法可以看做是二分法的扩展。若$f\left(t\right)$在区间$\left[a,b\right]$上连续且f$\left(a\right)\cdot f\left(b\right)<0$，则$[\left[a,b\right]$是含根区间。此时做$f(t)$过A$(a,f(a))$、$B\left(b,f\left(b\right)\right)$点的弦AB，则根据两点式有
$$L:y=\frac{b-x}{b-a}f(a)+\frac{x-a}{b-a}f(b)$$
弦AB与x轴交于
$$x=a-\frac{f(a)}{f(b)-f(a)}(b-a)$$
并作为$f\left(t\right)$在$\left[a,b\right]$上的根。不断重复上述过程直至含根区间足够小后停止迭代，与二分法不同的是，弦截法每次获得的新区间可能不是等分的，且通过固定b点后不断做弦，此时的搜索速度相较于二分法更快。弦截法的过程如下图所示

可见弦截法所得直线L与横轴交点总位于根的左侧，使得$f\left(x_n\right)\cdot f\left(x_0\right)<0$总是能够被满足。

x = np.arange(-1, 5, 0.01)
k = 0.7

def f(x):
    return np.exp(k*x)-6
def d_f(x):
    return k * exp(k * x)


x0 = 5.0; f0 = f(x0)
x1 = 0.3; f1 = f(x1)
n = 7

plt.figure(figsize=(10, 8))
plt.plot(x, f(x), linewidth=2.5, label="f(x)")
plt.vlines(x=0, ymin=f(min(x))-1, ymax=f(max(x))+1, color="black")
plt.hlines(y=0, xmin=-1, xmax=x0+0.5, color="black")
plt.vlines(x=x0, ymin=0, ymax=f(x0), linestyle="--")


marker_fx = []
marker_x = []
for i in range(1, n):
    lab = "$f(x_" + str(i) + ")$"
    lab2 = "$x_" + str(i) + "$"
    marker_fx.append(lab)
    marker_x.append(lab2)

for i in range(1, n):
    cut = x1 - (x0-x1)*(f(x1))/(f(x0)-f(x1))
    xx = np.arange(x1, 5, 0.01)
    yy = f(x1)*(x0-xx)/(x0-x1) + f(x0)*(xx-x1)/(x0-x1)
    plt.plot(xx, yy, color="black", linestyle="--", linewidth=0.8)
    plt.scatter(x1, f(x1)-1.3, c='black', marker=marker_fx[i-1], s=700 , edgecolors="black", alpha = 1, linewidths = 1) #fx
    plt.scatter(x1, 1.0, c='black', marker = marker_x[i-1], s=size_x , edgecolors="black", alpha = 1, linewidths = 0.5)   #x
    plt.vlines(x=x1, ymin=0, ymax=f(x1), linestyle="--")
    #x0 = x1
    x1 = cut
    
plt.yticks(fontstyle = "normal", fontproperties = 'Times New Roman', size = 13)
plt.xticks(fontstyle = "normal", fontproperties = 'Times New Roman', size = 13)
plt.scatter(x0, f(x0), c='yellow', marker = 'o', s=size_o , edgecolors="black", alpha = 1, linewidths = 1.5)
plt.scatter(x0+0.5, 0, c='black', marker = '>', s=100, edgecolors="black", alpha = 1, linewidths = 0.5)
plt.scatter(0, f(max(x))+1, c='black', marker = '^', s=100, edgecolors="black", alpha = 1, linewidths = 0.5)
plt.scatter(-0.2, f(max(x))+1, c='black', marker = '$f(x)$', s=600, edgecolors="black", alpha = 1, linewidths = 0.5)
plt.scatter(x0+0.5, -1.2, c='black', marker = '$x$', s=130, edgecolors="black", alpha = 1, linewidths = 0.5)
plt.scatter(4.85, 1, c='black', marker = '$x_0$', s=230, edgecolors="black", alpha = 1, linewidths = 0.5)
plt.legend(prop=fontdict_ita)
plt.tight_layout()
plt.show()

不动点法

不动点迭代法又称简单迭代。我们常用的牛顿法就是不动点迭代算法的特殊形式。为了求解方程$f\left(x\right)=0$，首先将$f\left(x\right)$转换为$x=g\left(x\right)$且需要满足$\left|g\left(x\right)\right|<1$，给定迭代初值$x_0$，循环计算$x_{n+1}=g\left(x_n\right)$直到$x_{n+1}-x_n$小于设定误差限$\epsilon$。然而该方法并不保证一定收敛至根附近，相反如果初值距离根过远还会导致结果发散。例如当我们想要求解方程$2x^3-x-1=0$的根，方程$f\left(x\right)=0$的转换方式不唯一，可以将其转化为$x=2x^3-1$或$x={\left(\left(x+1\right)/2\right)}^{1/3}$。从几何角度来看求解根的问题转变成为了求解$y=x$与$g\left(x\right)$交点问题。不断计算迭代式后就能得到根的近似值。下图展示不动点迭代算法的过程。

import numpy as np
x = np.arange(-5, 5, 0.01)
x0 = 2
def f(x):
    return 2*x**3-x-1

def g(x):
    return pow((x+1)/2,1/3)
def f2(x):
    return x
xk = []
for k in range(20):
    xk1 = g(x0)
    xk.append(g(x0))
    x0 = xk1
import matplotlib.pyplot as plt
fontdict_nor = {'family' : 'Times New Roman', 'size': 18, "style":"normal"}
fontdict_ita = {'family' : 'Times New Roman', 'size': 18, "style":"italic"}
plt.style.use("seaborn")
plt.style.use("seaborn")
x = np.arange(-0.5, 1.5, 0.01)

plt.figure(figsize=(10, 8))

n = len(xk)
marker_fx = []
marker_x = []
for i in range(1, n):
    lab = "$f(x_" + str(i) + ")$"
    lab2 = "$x_" + str(i) + "$"
    marker_fx.append(lab)
    marker_x.append(lab2)


plt.plot(x, f(x), label="f (x)", linewidth=2)
plt.plot(x, g(x), label="g (x)", linestyle="--")
plt.plot(x, f2(x),label="y=x", color="black")

plt.vlines(x=0, ymin=f(min(x))-1, ymax=f(max(x))+1, color="black")
plt.hlines(y=0, xmin=-0.5, xmax=1.5, color="black")

for i in range(len(xk)):
    plt.scatter(xk[i], 0, c='yellow', marker = 'o', s=80, edgecolors="black", alpha = 1, linewidths = 0.5)
    
for i in range(2):
    plt.scatter(xk[i], -0.2, marker = marker_x[i], s=200, edgecolors="black", alpha = 1, linewidths = 0.5)
    
plt.yticks(fontstyle = "normal", fontproperties = 'Times New Roman', size = 13)
plt.xticks(fontstyle = "normal", fontproperties = 'Times New Roman', size = 13)

plt.scatter(x0+0.5, 0, c='black', marker = '>', s=100, edgecolors="black", alpha = 1, linewidths = 0.5)
plt.scatter(0, f(max(x))+1, c='black', marker = '^', s=100, edgecolors="black", alpha = 1, linewidths = 0.5)
plt.scatter(-0.07, f(max(x))+1, c='black', marker = '$f(x)$', s=600, edgecolors="black", alpha = 1, linewidths = 0.5)
plt.scatter(1.5, -0.3, c='black', marker = '$x$', s=130, edgecolors="black", alpha = 1, linewidths = 0.5)
plt.scatter(xk[-1], f2(xk[-1]), c='yellow', marker = 'o', s=130, edgecolors="black", alpha = 1, linewidths = 0.5)


plt.vlines(x=xk[-1], ymin=0, ymax=f2(xk[-1]), color="black", linestyle="--")

plt.legend(prop=fontdict_ita)
plt.tight_layout()
plt.show()