牛顿—拉夫森迭代、加速牛顿—拉夫森迭代和哈利法改进牛顿——拉夫森迭代

天下弈星~

已于 2023-03-26 15:33:40 修改

阅读量1.2k

点赞数 3

分类专栏：数值方法文章标签：算法 python matlab 线性代数

于 2023-03-26 15:21:24 首次发布

本文链接：https://blog.csdn.net/m0_70452407/article/details/129778933

版权

数值方法专栏收录该内容

9 篇文章 0 订阅

订阅专栏

2.2加速牛顿—拉夫森迭代的matlab实现

2.2加速牛顿—拉夫森迭代的python实现

三、哈利法改进牛顿—拉夫森迭代

3.1哈利法改进牛顿—拉夫森迭代的原理

3.2哈利法改进的牛顿—拉夫森迭代的matlab实现

3.3哈利法改进的牛顿—拉夫森迭代的python实现

一、牛顿—拉夫森迭代

1.1牛顿—拉夫森定理

设 $f\in C^{2}[a,b]$ ，且存在数 $p\in [a,b]$ ，满足 $f(p)=0$ 。如果 $f^{'}(p)\neq 0$ ，则存在一个数 $\delta> 0$ ，对任意初始近似值 $p_{0}\in [p-\delta ,p+\delta ]$ ，使得由如下迭代定义的序列 $\left \{ p_{k}\right \}_{k=0}^{\bowtie }$ 收敛到p：

$p_{k}=g(p_{k-1})=p_{k-1}-\frac{f(p_{k-1})}{f^{'}(p_{k-1})},k=1,2,...$

其中，函数g(x)由如下公式定义：

$g(x)=x-\frac{f(x)}{f^{'}(x)}$

并被称为牛顿—拉夫森迭代公式。

1.2牛顿—拉夫森迭代的matlab实现

function [k,p0,err,ferr,y]=newton(f,df,p0,p,delta,epsilon,max1)
%输入  -f代表要求解的函数
%      -df代表所要求解的函数的一阶导数
%      -p0代表初始值
%      -p代表函数f的根
%      -M代表函数f的阶数
%      -delta代表p0所允许的误差
%      -epsilon代表函数值y所允许的误差
%      -max1代表迭代的次数
%输出  -k代表迭代的次数
%      -p0代表利用牛顿-拉夫森算法得到的近似值
%      -err代表利用牛顿-拉夫森算法得到的近似值与前一个近似值的差的绝对值
%      -ferr代表利用牛顿—拉夫森算法得到的近似值与根的差值
%      -y代表利用牛顿-拉夫森算法得到的近似值的函数值
for k=0:max1 
       p1=p0-feval(f,p0)/feval(df,p0);
       err=abs(p1-p0);
       relerr=2*err/(abs(p1)+delta);
       ferr=p-p0;
       y=feval(f,p0);
       disp([k,p0,err,ferr,y]);
       p0=p1;
       if (err<delta)||(relerr<delta)||(abs(y)<epsilon),break,end
end

在命令行窗口中输入：

>> format long

>> f=@(x) (x-1)*log(x);

>> df=@(x) log(x)+(x-1)/x;

>> newton(f,df,2,1,0,0,20)

最后可以得到结果如下：

k	P0	err	ferr	y
0	2.000000000000000	0.580940215803595	-1.000000000000000	0.693147180559945
1	1.419059784196405	0.227286602334084	-0.419059784196405	0.146668631586085
2	1.191773181862321	0.100028181255526	-0.191773181862321	0.033645121749720
3	1.091745000606795	0.046871620829905	-0.091745000606795	0.008053131572633
4	1.044873379776890	0.022681998076316	-0.044873379776890	0.001969748883255
5	1.022191381700575	0.011156462980106	-0.022191381700575	0.000487072782596
6	1.011034918720469	0.005532583287656	-0.011034918720469	0.000121102475330
7	1.005502335432813	0.002754940076509	-0.005502335432813	0.000030192705984
8	1.002747395356304	0.001374639691487	-0.002747395356304	0.000007537831277
9	1.001372755664817	0.000686613201220	-0.001372755664817	0.000001883165848
10	1.000686142463597	0.000343130057190	-0.000686142463597	0.000000470630039
11	1.000343012406407	0.000171520907450	-0.000343012406407	0.000000117637337
12	1.000171491498957	0.000085749425277	-0.000171491498957	0.000000029406813
13	1.000085742073679	0.000042871955757	-0.000085742073679	0.000000007351388
14	1.000042870117923	0.000021435288686	-0.000042870117923	0.000000001837808
15	1.000021434829236	0.000010717472049	-0.000021434829236	0.000000000459447
16	1.000010717357187	0.000005358692951	-0.000010717357187	0.000000000114861
17	1.000005358664236	0.000002679335707	-0.000005358664236	0.000000000028715
18	1.000002679328529	0.000001339665162	-0.000002679328529	0.000000000007179
19	1.000001339663367	0.000000669831908	-0.000001339663367	0.000000000001795
20	1.000000669831459	0.000000334915786	-0.000000669831459	0.000000000000449

从上面的结果可以明显地看出，牛顿—拉夫森迭代法的收敛速度很慢。

1.3牛顿—拉夫森迭代的python实现

import pandas as pd
import numpy as np
import math

# 方程式
def f (x) :
  return (x-1)*(math.log(x))

# 方程式的导数
def d_f (x) :
  return (math.log(x))+(x-1)/x

# 自行设定起始点x0
x = 2

# 迭代实现求解
d = {"x" : [x], "f(x)": [f(x)]}
for i in range(0, 20):
  x = x - f(x) / d_f(x)
  d["x"].append(x)
  d["f(x)"].append(f(x))

# 调用pandas库打印出x以及f(x)
data=pd.DataFrame(d, columns=['x', 'f(x)'])

print(data)

二、加速—牛顿拉夫森迭代

2.1牛顿—拉夫森迭代的加速收敛

设牛顿—拉夫森算法产生的序列线性收敛到M阶根x=p，其中M>1，则牛顿—拉夫森迭代公式 $p_{k}=p_{k-1}-\frac{Mf(p_{k-1})}{f^{'}(p_{k-1})}$

将产生一个收敛序列 $\left \{ p_{k}\right \}_{k=0}^{\bowtie }$ 二次收敛到p。

2.2加速牛顿—拉夫森迭代的matlab实现

function [k,p0,err,ferr,y]=newton(f,df,p0,p,M,delta,epsilon,max1)
%输入  -f代表要求解的函数
%      -df代表所要求解的函数的一阶导数
%      -p0代表初始值
%      -p代表函数f的根
%      -M代表函数f的阶数
%      -delta代表p0所允许的误差
%      -epsilon代表函数值y所允许的误差
%      -max1代表迭代的次数
%输出  -k代表迭代的次数
%      -p0代表利用加速牛顿-拉夫森算法得到的近似值
%      -err代表利用加速牛顿-拉夫森算法得到的近似值与前一个近似值的差的绝对值
%      -ferr代表利用加速牛顿—拉夫森算法得到的近似值与根的差值
%      -y代表利用加速牛顿-拉夫森算法得到的近似值的函数值
for k=0:max1 
       p1=p0-M*feval(f,p0)/feval(df,p0);
       err=abs(p1-p0);
       relerr=2*err/(abs(p1)+delta);
       ferr=p-p0;
       y=feval(f,p0);
       disp([k,p0,err,ferr,y]);
       p0=p1;
       if (err<delta)||(relerr<delta)||(abs(y)<epsilon),break,end
end

在命令行窗口中输入

>> format long

>> f=@(x) (x-1)*log(x);

>> df=@(x) log(x)+(x-1)/x;

>> newton(f,df,2,1,2,0,0,20)

最后得到的结果是：

k	P0	err	ferr	y
0	2.000000000000000	1.161880431607190	-1.000000000000000	0.693147180559945
1	0.838119568392810	0.154633333804432	0.161880431607190	0.028587194791166
2	0.992752902197242	0.007233911914516	0.007247097802758	0.000052711661389
3	0.999986814111759	0.000013185844774	0.000013185888241	0.000000000173869
4	0.999999999956533	0.000000000043467	0.000000000043467	0.000000000000000
5	1.000000000000000	0.000000000000000	0.000000000000000	0.000000000000000

很明显可以看出，加速牛顿—拉夫森算法迭代的速度比牛顿—拉夫森算法的迭代速度快很多。

2.2加速牛顿—拉夫森迭代的python实现

import pandas as pd
import numpy as np
import math

# 方程式
def f (x) :
  return (x-1)*(math.log(x))

# 方程式的导数
def d_f (x) :
  return (math.log(x))+(x-1)/x

# 自行设定起始点x0
x = 2

#确定阶数
M=2

# 迭代实现求解
d = {"x" : [x], "f(x)": [f(x)]}
for i in range(0, 5):
  x = x - M*f(x) / d_f(x)
  d["x"].append(x)
  d["f(x)"].append(f(x))

# 调用pandas库打印出x以及f(x)
data=pd.DataFrame(d, columns=['x', 'f(x)'])

print(data)

三、哈利法改进牛顿—拉夫森迭代

3.1哈利法改进牛顿—拉夫森迭代的原理

哈利（Hally）法是加速牛顿法收敛的另一个途径。哈利迭代公式是

$g(x)=x-\frac{f(x)}{f^{'}(x)}(1-\frac{f(x)f^{''}(x)}{2(f^{'}(x))^{2}})^{-1}$

括号中的项是对牛顿—拉夫森公式的改进。哈利法在f(x)的单根情况下可达到三次收敛（R=3）。

3.2哈利法改进的牛顿—拉夫森迭代的matlab实现

下面的matlab代码和python均是求解用其求解函数 $f(x)=x^{3}-3x+2$ 的单根， $p_{0}=-2.4$ 。

function [k,p0,err,y]=halley(f,df1,df2,p0,delta,epsilon,max1)
%输入  -f代表要求解的函数
%      -df1代表所要求解的函数的一阶导数
%      -df2代表所要求解的函数的二阶导数
%      -p0代表初始值
%      -delta代表p0所允许的误差
%      -epsilon代表函数值y所允许的误差
%      -max1代表迭代的次数
%输出  -k代表迭代的次数
%      -p0代表利用哈利法得到的近似值
%      -err代表利用哈利法得到的近似值与前一个近似值的差的绝对值
%      -y代表利用哈利法得到的近似值的函数值
for k=0:max1 
       p1=p0-(feval(f,p0)/feval(df1,p0))/(1-(feval(f,p0)*feval(df2,p0))/(2*feval(df1,p0)*feval(df1,p0)));
       err=abs(p1-p0);
       relerr=2*err/(abs(p1)+delta);
       y=feval(f,p0);
       disp([k,p0,err,y]);
       p0=p1;
       if (err<delta)||(relerr<delta)||(abs(y)<epsilon),break,end
end

接着，在命令行窗口中输入

>> format long;

>> f=@(x) x^3-3*x+2;

>> df1=@(x) 3*x^2-3;

>> df2=@(x) 6*x;

>> halley(f,df1,df2,-2.4,0.0000000001,0.0000000001,5)

最后可以得到以下结果：

k	P0	err	y
0	-2.400000000000000	0.386991869918699	-4.623999999999999
1	-2.013008130081301	0.013007408961696	-0.118090640545510
2	-2.000000721119605	0.000000721119605	-0.000006490079567
3	-2.000000000000000	0.000000000000000	0.000000000000000

3.3哈利法改进的牛顿—拉夫森迭代的python实现

import pandas as pd

# 方程式
def f (x) :
  return x**3-3*x+2

# 方程式的一阶导数
def d1_f (x) :
  return 3*x**2-3

#方程式的二阶导数
def d2_f(x):
    return 6*x

# 自行设定起始点x0
x = -2.4

# 迭代实现求解
d = {"x" : [x], "f(x)": [f(x)]}
for i in range(0, 3):
  x=x-f(x)/(d1_f(x)*(1-(f(x)*d2_f(x))/(2*d1_f(x)*d1_f(x))))
  d["x"].append(x)
  d["f(x)"].append(f(x))

# 调用pandas库打印出x以及f(x)
data=pd.DataFrame(d, columns=['x', 'f(x)'])

print(data)