一. 概述
五次及以上多项式方程没有根式解（就是没有像二次方程那样的万能公式），这个是被伽罗瓦用群论做出的最著名的结论。 那么这样的方程该如何求根呢???
牛顿于 1736 年公开提出了 牛顿迭代法来解决这个问题.
其核心思想 是 “逼近”.
二. 几何意义
我们知道，求一个方程的根，等同于求一个方程与x轴的交点的横坐标.
如何求交点的横坐标呢???
根据上图我们可以观察到:
从 x0 到 x1 到 x2, 它们离交点的横坐标越来越近了!!!
其最终会无限接近 交点的横坐标!!
迭代过程:
- 任选一点 x 0 x_0 x0 作垂线与 f(x) 交于 ( x 0 , f ( x 0 ) ) (x_0, f(x_0)) (x0,f(x0))
- 过 ( x 0 , f ( x 0 ) ) (x_0, f(x_0)) (x0,f(x0)) 作 f(x)图像的切线, 与 x 轴交于 x 1 x_1 x1
- 以 x1 为新点做垂线… 重复 1,2,3 步骤,直到得到满意的近似解 x n x_n xn
以上就是牛顿迭代法的全部过程了.
推论:
=> 斜率 =
f
′
(
x
0
)
=
f
(
x
0
)
x
0
−
x
1
f'(x0) = \frac{f(x0)}{x0-x1}
f′(x0)=x0−x1f(x0)
=>
x
1
=
x
0
−
f
(
x
0
)
f
′
(
x
0
)
x1 = x0 - \frac{f(x0)}{f'(x0)}
x1=x0−f′(x0)f(x0)
=>
x
n
+
1
=
x
n
−
f
(
x
n
)
f
′
(
x
n
)
x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}
xn+1=xn−f′(xn)f(xn)
三. 应用
3.1 求函数的根
def Dfunc(function, guess_value, delta=1e-5):
'''
通过导数的定义式求f(guess_value)这点的导数
'''
return (function(guess_value+delta)-function(guess_value))/delta
def newton_update(function, guess_value):
'''
通过迭代公式算出新的值
'''
return guess_value - function(guess_value)/Dfunc(function,guess_value)
def is_close_root(function,guess_value):
'''
判断是否已接近我们要求的根
'''
return function(guess_value) == 0 #精度要求过高可能会造成迭代次数过多
#return function(guess_value) < 1e-8 #折中方案
def newton_func(function, guess_value=2):
while not is_close_root(function,guess_value):
guess_value = newton_update(function, guess_value)
return guess_value
def find_root(function):
'''
lambda 作为函数传入，求函数的根
'''
return newton_func(function)
print(find_root(lambda x:pow(x,2)-9))#传入的是函数
print(find_root(lambda x:pow(x,2)-16))
print(find_root(lambda x:pow(x,3)-27))
print(find_root(lambda x:pow(x,3)-8))
3.2 开ｎ次方： k a k\sqrt{a} ka
要求 1 2 256 \frac{1}{2}\sqrt{256} 21256 即解方程 x 2 = 256 x^2=256 x2=256, 即求 x 2 − 256 = 0 x^2-256=0 x2−256=0 的根(移项即可)．
推广为: k a k\sqrt{a} ka ------> 解 x 1 k = a x^{\frac{1}{k}}=a xk1=a ------> 求 y = x 1 k − ａ 的 根 y=x^{\frac{1}{k}}-ａ的根 y=xk1−ａ的根
#求根只需要调用之前写的 求根函数即可
print(find_root(lambda x:pow(x,2)-256))
print(find_root(lambda x:pow(x,3)-27))
print(find_root(lambda x:pow(x,2)-2))
3.21 优化: 牛顿迭代方程的 化简
迭代方程为:
=>
x
n
+
1
=
x
n
−
f
(
x
n
)
f
′
(
x
n
)
x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}
xn+1=xn−f′(xn)f(xn)
当
f
(
x
n
)
=
x
2
−
a
时
f(x_n) = x^2 -a 时
f(xn)=x2−a时
=>
x
n
+
1
=
x
n
−
x
n
2
+
a
2
x
n
x_{n+1} = x_n - \frac{x_n^2+a}{2x_n}
xn+1=xn−2xnxn2+a
=>
x
n
+
1
=
2
x
n
2
−
x
n
2
+
a
2
x
n
x_{n+1} = \frac{2x_n^2-x_n^2 +a}{2x_n}
xn+1=2xn2xn2−xn2+a
=>
x
n
+
1
=
x
n
2
+
a
2
x
n
x_{n+1} = \frac{x_n^2+a}{2x_n}
xn+1=2xnxn2+a
=>
x
n
+
1
=
(
x
n
2
+
a
x
n
)
/
2
x_{n+1} = (\frac{x_n^2+a}{x_n})/2
xn+1=(xnxn2+a)/2
=>
x
n
+
1
=
(
x
n
+
a
x
n
)
/
2
x_{n+1} = (x_n + \frac{a}{x_n})/2
xn+1=(xn+xna)/2
对
f
(
x
)
=
x
n
−
a
f(x) = x^n -a
f(x)=xn−a 这样的函数的牛顿迭代求根公式可以化简为:
x
n
+
1
=
(
n
−
1
)
x
n
+
a
p
o
w
(
x
n
,
n
−
1
)
n
x_{n+1} = \frac{(n-1)x_n+\frac{a}{pow(x_n, n-1)}}{n}
xn+1=n(n−1)xn+pow(xn,n−1)a
化简后的迭代公式,更快!
#使用优化后迭代公式的代码,仅仅针对于 f(x) = x^n -a
def newton_update(base, expoent, guess_value):
return ((expoent-1)*guess_value + base/pow(guess_value, expoent-1))/expoent
def is_close_root(base, expoent, guess_value):
return abs(pow(guess_value, expoent) - base ) < 1e-12
def newton_func(base, expoent, guess_value):
while not is_close_root(base, expoent, guess_value):
guess_value = newton_update(base, expoent, guess_value)
return guess_value
def find_x_n_a(base, expoent, guess_value=2):
'''
求 y = x^n -a 的根
'''
return newton_func(base, expoent, guess_value)
def sqrt(base):
return find_x_n_a(base, 2, (base+1)/2)
print(find_x_n_a(256, 2))
print(find_x_n_a(2, 2))
print(find_x_n_a(8, 2))
print(find_x_n_a(8, 3))
cpp:
#include <iostream>
#include <cmath>
using namespace std;
double find_x_n_a(double base, double expoent, double guess_value=2){
return abs(pow(guess_value, expoent)-base)<1e-10 ? guess_value
: find_x_n_a(base, expoent,
((expoent-1)*guess_value + base/pow(guess_value, expoent-1))/expoent);
}
int main(){
cout << find_x_n_a(9,2) << endl;
cout << find_x_n_a(27,3) << endl;
cout << find_x_n_a(8,2) << endl;
cout << find_x_n_a(2,2) << endl;
return 0;
}
3.3 函数pow(base, expoent)
base^expoent 即 y = a x y = a^x y=ax 这个函数的求值操作
x为正数:
x 为整数，可以用 快速幂算法 求解
x 为小/分数，实际上是 开方操作
x 为负数:
对正数的结果 求倒数 即可. 或者： 求 ( 1 b a s e ) − x 求(\frac{1}{base})^{-x} 求(base1)−x 次方(更快)
开方 幂函数 指数函数 1 2 a \frac{1}{2}\sqrt{a} 21a a 1 2 a^{\frac{1}{2}} a21 y = x 2 − a 求 根 y = x^2-a求根 y=x2−a求根 1 3 a \frac{1}{3}\sqrt{a} 31a a 1 3 a^{\frac{1}{3}} a31 y = x 3 − a 求 根 y = x^3 -a 求根 y=x3−a求根 2 3 a \frac{2}{3}\sqrt{a} 32a a 2 3 a^{\frac{2}{3}} a32 y = x 3 2 − a 求 根 y = x^{\frac{3}{2}} -a 求根 y=x23−a求根 1 2 3 a \frac{1}{\frac{2}{3}\sqrt{a}} 32a1 a − 2 3 a^{-\frac{2}{3}} a−32 y = x 3 2 − a 的 正 根 的 倒 数 y = x^\frac{3}{2} - a 的正根的倒数 y=x23−a的正根的倒数
- 选 a 2 \frac{a}{2} 2a 为起始点, 其垂线交于 ( a 2 , f ( a 2 ) ) (\frac{a}{2}, f(\frac{a}{2})) (2a,f(2a))
- 求此点处的导数
- 用之前的推论算出起始点的下一个点 x n e x t x_{next} xnext
- 迭代多次得出近似解 x n x_n xn
写成 python 版本即为:
def newton_update(base, expoent, guess_value):
return ((expoent-1)*guess_value + base/pow(guess_value, expoent-1))/expoent
def is_close_root(base, expoent, guess_value):
return abs(pow(guess_value, expoent) - base ) < 1e-13
def newton_func(base, expoent, guess_value):
while not is_close_root(base, expoent, guess_value):
guess_value = newton_update(base, expoent, guess_value)
return guess_value
def find_x_n_a(base, expoent, guess_value=2):
'''
求 y = x^n -a 的根
'''
return newton_func(base, expoent, guess_value)
def int_quick_pow(base, expoent):
result = 1
while expoent:
if expoent & 1:
result *= base
expoent >>= 1
base *= base
return result
def pow_h(base, expoent):
'''
y = a^x y = base^expoent
'''
if expoent < 0:
'''
a^-p = (1/a)^p 或 y = x^(1/expoent)-a
'''
return pow_h(1/base, -expoent)# a(-2/3) = 1/(a^(2/3)) = (1/a)^(2/3)
elif expoent == 0:
return base
elif not expoent %1:
return int_quick_pow(base, expoent)
return find_x_n_a(base, 1/expoent, (base+1)/1) # y = x^(1/expoent) - a
'''
1 收敛缓慢或者不收敛 用 牛顿迭代法
'''
'''
小数倍 => 1.用泰勒级数
2.如果在 cpp 等静态类型语言中，由于 (base, fenzi, femu), n/m根号a 可以使用 m根号a^n 来计算
'''
print(pow_h(3,0))
print(pow_h(3,-1))
print(pow_h(4,-2))
print(pow_h(8,1/3))
print(pow_h(8,-1/3))
3.32 泰勒级数优化
快速收敛的幂级数要快于牛顿法
未完待续
6万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
