牛顿法的简单介绍及Matlab实现_matlab实现newton法

牛顿法原理简介

牛顿法的原理是利用函数$f(x)$ 的泰勒级数的前几项来寻找方程$f(x)=0$ 的根。

$f(x)$在$x_0$处的一阶泰勒展开
$$f(x)=f(x_0)+f^{\prime }(x_0)(x-x_0).$$
它的解
$$x_1=x_0-\frac{f(x_0)}{f^{\prime }(x_0)}.$$
进而推出
$$x_{n+1}=x_n-\frac{f(x_n)}{f^{\prime }(x_n)}.$$
由此迭代下去求出近似解。

参考百度百科牛顿迭代法.

使用牛顿法求解一元方程

使用牛顿法来解方程
$$2=(x-1{)}^2+2.$$

使用牛顿法求解，先给一个初始值$x_{guess}=10$，迭代结果不断逼近真实值。

迭代50次的结果与真实值，可以看出，牛顿法迭代逼近真实值的速度非常快，但并不是迭代次数越多越好。

代码如下

%牛顿法求一元二次方程的根 by wuyc
%y=(x-1)^2+2
%y'=2x-2

clc; clear all;
xguess = 10;%猜的初始点
N=50;%迭代次数

x1=xguess;
x=zeros(N,1);

for i=1:N
deltax= -((x1-1)^2+2)/(2*x1-2);%delta x=-f'(xn)/f(xn)
x2 = x1+deltax; %xn+1=xn+delta x
x1=x2;
x(i)=x1;%x是每次迭代的结果
end

figure (1);

t=[0:0.1:12];
y=(t-1).^2+2;

yx=(x-1).^2+2;

plot(t,y,'LineWidth',2);hold on; %画方程曲线

plot([xguess,xguess],[0,(xguess-1).^2+2]);
plot([xguess,x(1)],[(xguess-1).^2+2,0]);

for ii=1:2 
plot([x(ii),x(ii)],[0,yx(ii)]);
plot([x(ii),x(ii+1)],[yx(ii),0]);
end
set(gca,'LineWidth',2);
xlabel('x','FontSize',18);
ylabel('y','FontSize',18);

figure (2);

plot (x);
hold on
plot ([1,50],[1,1]);%真实值
set(gca,'LineWidth',2);
xlabel('迭代次数');
ylabel('迭代结果');
legend('牛顿法结果','真实值');
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47

使用牛顿法求解平面定位问题

定位问题应用十分广泛，这里我们把问题简化为二维的矩形区域定位。
由四个台站TZ1、TZ2、TZ3、TZ4决定的矩形区域内的某一点mt$(m_{true})$发出信号，信号的速度是$v$，那么根据四个台站接收到的时间的平方dobs($d_{observation})$，我们可以求出信号的位置。(台站只能记录信号出现的时刻，而不能记录到信号传播到台站所需的时间$t$，这里我们又把问题简化了。)
这里插一句，在反演问题中，我们常把模型参数称为$m$，数据称为$d$，所谓反演问题，就是给定数据$d$求解模型参数$m$的过程。
因为时间$t$等于路程$s$除以速度$v$，有
$$t^2=\frac{1}{v^2}\times s^2$$
因此对应的反演问题为
[ d 1 d 2 d 3 d 4 ] = 1 v 2 [ ( x 1 2 − m t 1 2 ) + ( y 1 2 − m t 2 2 ) ( x 2 2 − m t 1 2 ) + ( y 2 2 − m t 2 2 ) ( x 3 2 − m t 1 2 ) + ( y 3 2 − m t 2 2 ) ( x 4 2 − m t 1 2 ) + ( y 4 2 − m t 2 2 ) ]

$$\begin{bmatrix} d_1 \\ d_2 \\ d_3 \\ d_4 \\ \end{bmatrix}$$ =\frac{1}{v^2} $$\begin{bmatrix} (x_1^2-mt_1^2)+(y_1^2-mt_2^2) \\ (x_2^2-mt_1^2)+(y_2^2-mt_2^2) \\ (x_3^2-mt_1^2)+(y_3^2-mt_2^2) \\ (x_4^2-mt_1^2)+(y_4^2-mt_2^2) \\ \end{bmatrix}$$ ​d1​d2​d3​d4​​ ​=v21​ ​(x12​−mt12​)+(y12​−mt22​)(x22​−mt12​)+(y22​−mt22​)(x32​−mt12​)+(y32​−mt22​)(x42​−mt12​)+(y42​−mt22​)​ ​

这里的参数和问题显得有些混乱，主要目的是为了是反演问题简单。简单的说就是操场四个顶点有四个接收器，操场中间某一个点发出一个速度为$v$信号，发出信号的同时接收器开始计时，因此接收器记录的时刻$t$就是信号传输的时间。
我们要根据$t$来求解信号的位置$m$，为了使反演问题简单，我们把$t^2$看作数据$d$.

代码运行结果如下，黑点是猜的初始点，绿点是真实点，红点是反演结果。

代码如下

% 根据到时反演信号位置

clc;clear all;

N=4;
v=100;%速度
%4个台站的位置
TZ1=[0,0];
TZ2=[400,0];
TZ3=[400,600];
TZ4=[0,600];
TZ=[TZ1;TZ2;TZ3;TZ4];
x=TZ(:,1);%台站位置横坐标
y=TZ(:,2);%台站位置纵坐标

M=2;%反演的模型参数维度
mt = [100, 100]';%真实的信号位置

dtrue = (1/v^2).*[(x-ones(N,1)*mt(1)).^2 + (y-ones(N,1)*mt(2)).^2];%真实的到时的平方
sd=0.4;
dobs = dtrue + random('Normal',0,sd,N,1);%真实值 加上均值是0，方差是sd，N行1列的随机数(随机误差) 生成观测值（到时平方）

mg=[400,600]';%mguess 牛顿法需要一个初始点
dg = (1/v^2).*[(x-ones(N,1)*mg(1)).^2 + (y-ones(N,1)*mg(2)).^2];%dg是猜的点到时的平方
Eg = (dobs-dg)'*(dobs-dg);


%生成401*601网格
L = 401;
J = 601
Dm = 1;
m1min=0;
m2min=0;
ma1 = m1min+Dm*[0:L-1]';
ma2 = m2min+Dm*[0:J-1]';
m1max = ma1(L);
m2max = ma2(J);

% 计算误差
E = zeros(L,J);
for j = [1:L]
for k = [1:J]
    dpre = (1/v^2).*[(x-ones(N,1)*ma1(j)).^2 + (y-ones(N,1)*ma2(k)).^2];%预测值
    E(j,k) = (dobs-dpre)'*(dobs-dpre);%误差
end
end

figure(1);
%误差图
set(gca,'LineWidth',2);
hold on;
axis( [m2min, m2max, m1min, m1max] );
imagesc( [m2min, m2max], [m1min, m1max], E);
colorbar;
xlabel('m2','FontSize',18);
ylabel('m1','FontSize',18);
plot( mt(2), mt(1), 'go', 'LineWidth', 3 );%绿色点表示真实值
plot( mg(2), mg(1), 'ko', 'LineWidth', 3 );%黑色是猜的初始点

% Newton's method, calculate derivatives
% y = 1/v^2*((x-m1)^2+(y-m2)^2)
% dy/dm1 = 2/v^2*(m1-x)
% dy/dm2 = 2/v^2*(m2-y)

%Ehis用来储存每次迭代的误差 m1his、m2his用来储存每次迭代的结果
Niter=10;%迭代次数
Ehis=zeros(Niter+1,1);
m1his=zeros(Niter+1,1);
m2his=zeros(Niter+1,1);
Ehis(1)=Eg;
m1his(1)=mg(1);
m2his(1)=mg(2);

G = zeros(N,M);%G矩阵的维度是N*M
for k = [1:Niter]

    dg = (1/v^2).*[(x-ones(N,1)*mg(1)).^2 + (y-ones(N,1)*mg(2)).^2];%dg是猜的点到时的平方
    dd = dobs-dg;%猜的值与观测值的差
    Eg=dd'*dd;
    Ehis(k+1)=Eg;
    
    G = zeros(N,M);%导数矩阵
    G(:,1) = 2/v^2.*(ones(N,1)*mg(1)-x);
    G(:,2) = 2/v^2.*(ones(N,1)*mg(2)-y);
    
    
    dm = ((G'*G)^-1)*G'*dd;%最小二乘求解 dm 相当于delta x
    
    mg = mg+dm;
    plot( mg(2), mg(1), 'wo', 'LineWidth', 2 );%画出迭代后的点
    
    m1his(k+1)=mg(1);
    m2his(k+1)=mg(2);
    plot( [m2his(1+k-1), m2his(1+k) ], [m1his(1+k-1), m1his(1+k) ], 'r', 'LineWidth', 2 );
    
end


plot( m2his(Niter+1), m1his(Niter+1), 'ro', 'LineWidth', 2 );%红色是迭代最终结果
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99

此代码参考了William Menke 的Geophysical Data Analysis: Discrete Inverse Theory.

需要指出的是，虽然牛顿法原理简单收敛快，但这个过程并非总能实现，这里就涉及到全局极小和局部极小的问题了。幸运的是上述两个例子不存在这个问题，但在下图中，只有当初始值在红色区域内时，牛顿法才可以收敛到全局极小。

参考文献

链接：
1:牛顿迭代法.
2: Geophysical Data Analysis: Discrete Inverse Theory