牛顿法及其变形方法
牛顿法
牛顿法的构造
如果方程
f
(
x
)
=
0
f(x)=0
f(x)=0 是线性方程,则它的根是容易求解的. 牛顿法是一种线性化的近似方法,其基本思想是将非线性方程
f
(
x
)
=
0
f(x)=0
f(x)=0转化为线性方程来迭代求解.
设
f
(
x
)
f(x)
f(x)在
x
∗
x^*
x∗ 附近二次连续可微, 是附近的一个近似值(假定
f
′
(
x
0
)
≠
0
f'(x_0)\neq0
f′(x0)=0 ),将函数在点
x
0
x_0
x0处做一阶 Taylor 展开,有
f
(
x
)
≈
f
(
x
0
)
+
f
′
(
x
0
)
(
x
−
x
0
)
f(x)\approx f(x_0)+f'(x_0)(x-x_0)
f(x)≈f(x0)+f′(x0)(x−x0)则方程
f
(
x
)
=
0
f(x)=0
f(x)=0近似为如下线性方程
f
(
x
0
)
+
f
′
(
x
0
)
(
x
−
x
0
)
=
0
f(x_0)+f'(x_0)(x-x_0)=0
f(x0)+f′(x0)(x−x0)=0其根记为
x
1
=
x
0
−
f
(
x
0
)
f
′
(
x
0
)
x_1=x_0-\frac{f(x_0)}{f'(x_0)}
x1=x0−f′(x0)f(x0),即为
x
∗
x^*
x∗的新近似值. 重复上述过程,将函数
f
(
x
)
f(x)
f(x)在点
x
1
x_1
x1展开得
f
(
x
)
≈
f
(
x
1
)
+
f
′
(
x
1
)
(
x
−
x
1
)
f(x)\approx f(x_1)+f'(x_1)(x-x_1)
f(x)≈f(x1)+f′(x1)(x−x1)则方程
f
(
x
)
=
0
f(x)=0
f(x)=0 近似为线性方程
f
(
x
1
)
+
f
′
(
x
1
)
(
x
−
x
1
)
=
0
f(x_1)+f'(x_1)(x-x_1)=0
f(x1)+f′(x1)(x−x1)=0
得到新的近似解
x
2
=
x
1
−
f
(
x
1
)
f
′
(
x
1
)
x_2=x_1-\frac {f(x_1)}{f'(x_1)}
x2=x1−f′(x1)f(x1)
如此继续下去,得到迭代序列
x
k
+
1
=
x
k
−
f
(
x
k
)
f
′
(
x
k
)
,
k
=
0
,
1
,
2
,
⋯
x_{k+1}=x_k -\frac{f(x_k)}{f'(x_k)},k=0,1,2,\cdots
xk+1=xk−f′(xk)f(xk),k=0,1,2,⋯
这种方法称为牛顿法,上式称为牛顿迭代公式.
几何解释
牛顿法在几何上是运用曲线 y = f ( x ) y=f(x) y=f(x)的切线与 x x x轴的交点坐标来近似曲线与 x x x轴交点坐标的,因此牛顿法也称为切线法。
计算步骤
函数代码
import random
import sympy
def Newton_itera(a, b, f, R):
x = sympy.symbols('x')
X_k = [random.uniform(a, b), "x2"] #构造XK序列
X_k[-1] = float(X_k[-2] - f.subs(x, X_k[-2]) / (sympy.diff(f, x).subs(x, X_k[-2])))
while (abs(X_k[-1] - X_k[-2]) >= R):
x_k = float(X_k[-1] - f.subs(x, X_k[-1]) / (sympy.diff(f, x).subs(x, X_k[-1])))
X_k.append(x_k)
return X_k
例题
我们可以通过牛顿法求解分数和根式
- C \sqrt{C} C
- 1 C \frac{1}{C} C1
当
C
=
115
时,求解
115
的近似值
C=115时,求解\sqrt{115}的近似值
C=115时,求解115的近似值
解 由于
C
\sqrt{C}
C是方程
x
2
−
C
=
0
x^2-C=0
x2−C=0的正根,因此取
f
(
x
)
=
x
2
−
C
,
f
′
(
x
)
=
2
x
f(x)=x^2-C,f'(x)=2x
f(x)=x2−C,f′(x)=2x构造迭代公式
x
k
+
1
=
1
2
(
x
k
+
C
x
k
)
,
k
=
0
,
1
,
2
,
⋯
x_{k+1}=\frac{1}{2}(x_k+\frac{C}{x_k}),k=0,1,2,\cdots
xk+1=21(xk+xkC),k=0,1,2,⋯
计算结果如下
k | 0 | 1 | 2 | 3 | 4 |
---|---|---|---|---|---|
x k x_k xk | 10 | 10.750000 | 10.723837 | 10.723805 | 10.723805 |
def Newton_itera(a, b, f, R):
x = sympy.symbols('x')
X_k = [10, "x2"]
X_k[-1] = float(X_k[-2] - f.subs(x, X_k[-2]) / (sympy.diff(f, x).subs(x, X_k[-2])))
while (abs(X_k[-1] - X_k[-2]) >= R):
x_k = float(X_k[-1] - f.subs(x, X_k[-1]) / (sympy.diff(f, x).subs(x, X_k[-1])))
X_k.append(x_k)
return X_k
def main():
x = sympy.symbols('x')
f = x**2-115
print(Newton_itera(10,11,f,10**-7))
if __name__ == "__main__":
main()
结果
[10, 10.75, 10.723837209302326, 10.723805294811097, 10.723805294763608]
牛顿法的变形
简化牛顿法
为了避免牛顿迭代法中的导数计算,取某常数
M
≠
0
M\neq0
M=0代替牛顿迭代公式中的
f
′
(
x
k
)
f'(x_k)
f′(xk),得到简化牛顿迭代公式
x
k
+
1
=
f
(
x
k
)
−
x
k
M
,
k
=
0
,
1
,
2.
⋯
x_{k+1}=f(x_k)-\frac{x_k}{M},k=0,1,2.\cdots
xk+1=f(xk)−Mxk,k=0,1,2.⋯显然,
M
最好取与
f
′
(
x
∗
)
M最好取与f'(x^*)
M最好取与f′(x∗)较为接近的常数,有时
M
取
f
′
(
x
0
)
M取f'(x_0)
M取f′(x0),此时简化牛顿法为线性收敛。
简化牛顿法的几何解释:过曲线
y
=
f
(
x
)
y=f(x)
y=f(x)上点
P
k
(
x
k
,
f
(
x
k
)
)
P_k(x_k,f(x_k))
Pk(xk,f(xk))作斜率
M
M
M的平行弦,与
x
x
x轴的交点横坐标
x
k
+
1
x_{k+1}
xk+1作
x
∗
x^*
x∗的近似值,因此简化牛顿法也称为平行弦法
牛顿下山法
由牛顿法局部收敛性可知,初始值 x 0 x_0 x0应选取根在 x ∗ x^* x∗附近,但对于有些问题,很难检验满足条件的初始值 x 0 x_0 x0,此时可利用下山法扩大初始值的选取范围,将牛顿迭代公式修改为 x k + 1 = x k − λ f ( x k ) f ′ ( x k ) , k = 0 , 1 , 2 , ⋯ x_{k+1}=x_k-\lambda\frac{f(x_k)}{f'(x_k)},k=0,1,2,\cdots xk+1=xk−λf′(xk)f(xk),k=0,1,2,⋯式中 λ \lambda λ是一个参数,其选取应使得 ∣ f ( x k + 1 ) ∣ < ∣ f ( x k ) ∣ |f(x_{k+1})|<|f(x_k)| ∣f(xk+1)∣<∣f(xk)∣成立,这种方法称为牛顿下山法, λ \lambda λ称为下山因子,满足 0 < ε r ≤ λ ≤ 1 0<\varepsilon_r\leq\lambda\leq1 0<εr≤λ≤1, ε r \varepsilon_r εr为下山因子下界,为方便,一般开始时可简单的取 λ = 1 \lambda = 1 λ=1,然后逐步分半减小,即选取 λ = 1 , 1 2 . 1 4 , ⋯ , λ ≥ ε r \lambda=1,\frac{1}{2}.\frac{1}{4},\cdots,\lambda \geq \varepsilon_r λ=1,21.41,⋯,λ≥εr,直至 ∣ f ( x k + 1 ) ∣ < ∣ f ( x k ) ∣ |f(x_{k+1})|<|f(x_k)| ∣f(xk+1)∣<∣f(xk)∣成立。
割线法
为避免计算f"(x),由导数定义,有
f
′
(
x
)
≈
f
(
x
k
)
−
f
(
x
k
−
1
)
x
k
−
x
k
−
1
f'(x)\approx\frac{f(x_k)-f(x_{k-1})}{x_k-x_{k-1}}
f′(x)≈xk−xk−1f(xk)−f(xk−1),代入牛顿迭代公式得
x
k
+
1
=
x
k
−
f
(
x
k
)
f
(
x
k
)
−
f
(
x
k
−
1
)
(
x
k
−
x
k
−
1
)
x_{k+1}=x_k-\frac{f(x_k)}{f(x_k)-f(x_{k-1})}(x_k-x_{k-1})
xk+1=xk−f(xk)−f(xk−1)f(xk)(xk−xk−1)使用式(2.17)求解非线性方程的方法称为割线法(或弦截法).
割线法的几何解释:过曲线
y
=
f
(
x
)
y=f(x)
y=f(x)上点
P
k
−
1
(
x
k
−
1
,
f
(
x
k
−
1
)
)
P_{k-1}(x_{k-1},f(x_{k-1}))
Pk−1(xk−1,f(xk−1))及
P
k
(
x
k
,
f
(
x
k
)
)
P_k(x_k,f(x_k))
Pk(xk,f(xk))作割线,记割线与x轴交点的横坐标为
x
k
+
1
x_{k+1}
xk+1,用割线与
x
x
x轴的交点近似切线与x轴的交点.
利用割线法计算
x
k
+
1
x_{k+1}
xk+1,时要用到前两步的结果
x
k
x_k
xk,和
x
k
−
1
x_{k-1}
xk−1,这类算法称为两步迭代法.而牛顿法在计算
x
k
+
1
x_{k+1}
xk+1,时只用到前一步的值
x
k
x_k
xk,这类算法称为单步方法.因此,在使用割线法公式时必须先给出两个初始值
x
0
x_0
x0和
x
1
x_1
x1.
多项式方程求根法
很多实际应用问题需要求多项式的全部零点,它等价于求多项式方程( n>1 )
P
(
x
)
=
a
0
x
n
+
a
1
x
n
−
1
+
⋯
+
a
n
−
1
x
+
a
n
=
0
(2.20)
P(x)= a_0x^n+a_1x^{n-1}+\cdots+a_{n-1}x+a_n=0\tag{2.20}
P(x)=a0xn+a1xn−1+⋯+an−1x+an=0(2.20)
的全部根,式中系数
a
0
,
a
1
,
⋯
,
a
n
a_0,a_1,\cdots,a_n
a0,a1,⋯,an是实数.
前面介绍的任何非线性方程求根方法均适用于多项式方程求根,但由于多项式方程的特殊性,因此可以针对其特点提供更有效的算法,通常使用牛顿法最好.本节介绍求解多项式方程( n>1)的牛顿法.
定理:对于任意多项式
c
(
x
)
c(x)
c(x)和
d
(
x
)
d(x)
d(x),其中
d
(
x
)
≠
0
d(x)\neq0
d(x)=0,存在唯一的多项式
Q
(
x
)
Q(x)
Q(x)和
r
(
x
)
r(x)
r(x),满足
c
(
x
)
=
Q
(
x
)
d
(
x
)
+
r
(
x
)
c(x)=Q(x)d(x)+r(x)
c(x)=Q(x)d(x)+r(x)且
r
(
x
)
r(x)
r(x)的次数小于
d
(
x
)
d(x)
d(x)的次数.
Q
(
x
)
Q(x)
Q(x)称为商多项式,
r
(
x
)
r(x)
r(x)称为余多项式.
由定理可知,当
c
(
x
)
c(x)
c(x)为
n
n
n次,
d
(
x
)
d(x)
d(x)为
m
m
m次,且
m
<
n
m<n
m<n,
Q
(
x
)
Q(x)
Q(x)为
n
−
m
n-m
n−m次,
r
(
x
)
r(x)
r(x)次数小于
m
m
m次.
由于牛顿迭代公式中涉及函数和函数的导函数运算,因此,对于代数多项式方程,其计算量较大.多项式的牛顿法基本思想是:将牛顿法应用于求解多项式方程,再利用秦九韶方法求出牛顿公式中的函数值
f
(
x
k
)
f(x_k)
f(xk)和
f
′
(
x
k
)
f'(x_k)
f′(xk) .
将牛顿法应用于多项式方程,有迭代公式
x
k
+
1
=
x
k
−
P
n
(
x
k
)
P
n
′
(
x
k
)
,
k
=
=
0
,
1
,
2
,
⋯
(2.21)
x_{k+1}=x_k-\frac{P_n(x_k)}{P_n'(x_k)},k= =0,1,2,\cdots\tag{2.21}
xk+1=xk−Pn′(xk)Pn(xk),k==0,1,2,⋯(2.21)
借助秦九韶方法,不直接使用牛顿迭代公式以及其导数表达式,分别建立
P
(
x
k
)
P(x_k)
P(xk),
P
′
(
x
k
)
P'(x_k)
P′(xk)与多项式系数的关系,得到
P
n
(
x
k
)
P_n(x_k)
Pn(xk),
P
n
′
(
x
k
)
P_n'(x_k)
Pn′(xk),从而简化计算.
用
(
x
−
x
k
)
(x-x_k)
(x−xk)除多项式
P
n
(
x
)
P_n(x)
Pn(x),设商
Q
(
x
)
=
b
0
x
n
−
1
+
b
1
x
n
−
2
+
⋯
+
b
n
−
2
x
+
b
n
−
1
(2.22)
Q(x)=b_0x^{n-1}+b_1x^{n-2}+\cdots+b_{n-2}x+b_{n-1}\tag{2.22}
Q(x)=b0xn−1+b1xn−2+⋯+bn−2x+bn−1(2.22)
余数为
b
n
b_n
bn,,则
P
n
(
x
)
=
Q
(
x
)
(
x
−
x
k
)
+
b
n
(2.23)
P_n(x)=Q(x)(x -x_k)+b_n\tag{2.23}
Pn(x)=Q(x)(x−xk)+bn(2.23)
将多项式方程的基本形式和
Q
(
x
)
Q(x)
Q(x)带入上式得到如下递推关系
{
b
0
=
a
0
b
1
=
a
1
+
b
i
−
1
x
k
,
i
=
1
,
2
,
⋯
,
n
(2.24)
\begin{cases} b_0=a_0\\ b_1 = a_1+b_{i-1}x_k,&i = 1,2,\cdots,n \\ \end{cases}\tag{2.24}
{b0=a0b1=a1+bi−1xk,i=1,2,⋯,n(2.24)
同理,再用(x一x,)除Q(x),得到商
H
(
x
)
=
c
0
x
n
−
2
+
c
1
x
n
−
3
+
⋯
+
c
n
−
3
x
+
c
n
−
2
(2.25)
H(x)= c_0x^{n-2}+c_1x^{n-3}+\cdots+c_{n-3}x +c_{n-2}\tag{2.25}
H(x)=c0xn−2+c1xn−3+⋯+cn−3x+cn−2(2.25)
余数为
c
n
−
1
c_{n-1}
cn−1,则
Q
(
x
)
=
H
(
x
)
(
x
−
x
k
)
+
c
n
−
1
(2.26)
Q(x)= H(x)(x -x_k)+c_{n-1}\tag{2.26}
Q(x)=H(x)(x−xk)+cn−1(2.26)
将式(2.22)和式(2.25)代入上式,比较系数,得如下递推关系式
{
c
0
=
b
0
c
i
=
b
i
+
c
i
−
1
x
k
,
i
=
1
,
2
,
…
,
n
一
1
(2.27)
\begin{cases} c_0= b_0\\ c_i = b_i+c_{i-1}x_k , &i =1,2,…,n一1\\ \end{cases}\tag{2.27}
{c0=b0ci=bi+ci−1xk,i=1,2,…,n一1(2.27)
由式可得
Q
(
x
k
)
=
c
n
−
1
Q(x_k)=c_{n-1}
Q(xk)=cn−1,对基本形式求导,可得到
P
′
(
x
k
)
=
Q
(
x
k
)
=
c
n
−
1
P'(x_k)=Q(x_k)=c_{n-1}
P′(xk)=Q(xk)=cn−1因此,通过递推式(2.27)计算
c
n
−
1
c_{n-1}
cn−1,可得到
P
n
′
(
x
k
)
P'_n(x_k)
Pn′(xk).
计算步骤
综上所述,牛顿法求多项式方程根的计算步骤归纳如下.
- 取初始值 x 0 x_0 x0,通常可以取 x 0 = 0 x_0=0 x0=0.
- 对
k
=
0
,
1
,
2
,
…
k =0,1,2,…
k=0,1,2,…计算
b 0 = a 0 , b i = a i + b i − 1 x k , ( i = 1 , 2 , ⋯ , n ) b_0=a_0,b_i=a_i+b_{i-1}x_k ,(i=1,2,\cdots,n) b0=a0,bi=ai+bi−1xk,(i=1,2,⋯,n) c 0 = b 0 , c j = b j + c j − 1 x k , ( j = 2 , 3 , ⋯ , n − 1 ) c_0=b_0,c_j=b_j+c_{j-1}x_k,(j=2,3,\cdots,n-1) c0=b0,cj=bj+cj−1xk,(j=2,3,⋯,n−1) x k + 1 = x k − b n c n − 1 x_{k+1} = x_k - \frac{b_n}{c_{n-1}} xk+1=xk−cn−1bn - 当误差 ∣ b n ∣ < ε |b_n|<\varepsilon ∣bn∣<ε或 ∣ x k + 1 一 x k ∣ = ∣ b n c n − 1 ∣ < ε |x_{k+1}一x_k|=|\frac{b_n}{c_{n-1}}|<\varepsilon ∣xk+1一xk∣=∣cn−1bn∣<ε时,迭代终止;否则,转步骤2继续迭代.
代码实现
def multinomial_root(coef,R,x0):
'''多项式求根公式的迭代实现'''
#读入多项式阶数和迭代初值x0
n= len(coef)
X_k = [x0]
while True:
B, C = [coef[0]], [coef[0]] #初始化b0 = c0 =a0
for i in range(1, n):
B.append(coef[i] + B[i - 1] * X_k[-1]) #计算bi
C.append(B[i] + C[i - 1] * X_k[-1]) #计算cj
X_k.append(X_k[-1] - B[n - 1] / C[n - 2]) #计算x{k+1}
if (abs(X_k[-1]-X_k[-2]))< R: # 判断 |x{k+1}-xk| 是否达到精度
return X_k