【龙年开工：非线性方程（组）求解】牛顿-拉夫逊算法多输入单输出

最新推荐文章于 2024-10-31 16:16:13 发布

sanfan66

最新推荐文章于 2024-10-31 16:16:13 发布

阅读量1k

点赞数 20

分类专栏：优化理论及应用文章标签：算法学习

本文链接：https://blog.csdn.net/sanfan66/article/details/136140679

版权

优化理论及应用专栏收录该内容

5 篇文章 1 订阅

订阅专栏

文章目录

【非线性方程求解】割线法。

1. 概述

总结上篇博客内容，求解单变量非线性方程 $f (x) = 0$ 时，牛顿法和割线法的区别在于当前导数的求解方法不同，之后都是用一条直线去近似 $f (x)$ 。若求解的是非线性方程组，即含有多个输入和多个输出，实现有些不同。下面介绍牛顿法的具体实现，有时也叫作牛顿-拉夫逊（拉弗森）法。

2. 原理步骤

2.1 原理

输入 $\{x_{i}|i=1,\dots,nx\}$ 经过黑盒子后得到输出 $\{y_{j}|j=1,\dots,ny\}$

但当我想要某种结果 $\{y_{j}^{t}|j=1,\dots,ny\}$ 时，如何确定输入？为了解决这个问题，先按数学老师教的约定下符号：
$X=[x_{1},\dots,x_{nx}]^{T}$
$T$ 为转秩符号，表示 $X$ 是一个列向量。同理：
$\begin{aligned} Y=[y_{1},\dots,y_{ny}]^{T} \\ \\ Y^{t}=[y_{1}^{t},\dots,y_{ny}^{t}]^{T} \end{aligned}$
$t$ 表示目标输出。为了和之前讲过的步骤比较，把问题转换成求解以下方程，当然，你不想转换也行。
$f(X)=\varphi(X)-Y^{t}=0$
接下来使用以前的老套路，仅仅把 $x$ 换成 $X$ 。为避免与 $X$ 中的元素混淆，我把第 $k$ 次迭代后的结果写成上标：
$\begin{align} &\to f(X)=f\left(X^{k}\right)+f^{\prime}\left(X^{k}\right)\left(X-X^{k} \right)+{f^{\prime \prime}(\xi)}\left(X-X^{k}\right)^{2}/ (2!) \nonumber\\ \nonumber\\ &\to f\left(X^{k}\right)+f^{\prime}\left(X^{k}\right)\left(X-X^{k}\right)=0 \nonumber\\\nonumber\\ & \to f^{\prime}\left(X^{k}\right)\left(X-X^{k}\right)=-f\left(X^{k}\right) \\ \nonumber\\ &\to \cancel{X-X^{k}=-\frac{f\left(X^{k}\right)}{f^{\prime}\left(X^{k}\right)}} \nonumber\\ \nonumber\\ &\to \cancel{X=X^{k}-\frac{f\left(X^{k}\right)}{f^{\prime}\left(X^{k}\right)}} \nonumber \end{align}$

你不能继续套路下去了，因为带有 $X$ 的都是矩阵， $f^{\prime}\left(X^{k}\right)$ 表示 $f (x)$ 对 $X$ 的偏导数矩阵，需要根据 $f^{\prime}\left(X^{k}\right)$ 的维度作不同变换，怎么检查 $f^{\prime}\left(X^{k}\right)$ 写的对不对呢？暂且将导数理解为分母变化量为1时分子的变化量，则当 $X$ 变化量为 $d X$ 时 $f (X)$ 变化量为 $df (X)$ ，检查下式是否成立即可。

当系统为多输入单输出时，即 $n x > 1, n y = 1$ ，如果这是线性方程，直接用线性代数的方法求解公式（1）就行了，仅需一步并且能证明有很多解。对于非线性方程，什么情况都有可能存在，需要迭代。公式（1）展开为：
$\begin{equation} \frac{\partial f_1}{\partial x_1}\cdot dx_1+\frac{\partial f_1}{\partial x_2}\cdot dx_2+\cdots+\frac{\partial f_1}{\partial x_{nx}}\cdot dx_{nx}=df_1 \end{equation}$
我们看到，每个 $x_i$ 变化都会对 $f_1$ 产生影响，但问题我只想知道 $f_1=0$ 时 $X = ？$ 仿照公式（2）的形式，可以写出 $X$ 对 $f_1$ 的导数矩阵 $X^{\prime}_{f_1}$ 并且有：
$\begin{equation} \frac{\partial x_1}{\partial f_1}\cdot df_1+\frac{\partial x_2}{\partial f_1}\cdot df_1+\cdots+\frac{\partial x_{nx}}{\partial f_1}\cdot df_1=dX \end{equation}$
公式(3)表示 $f_1$ 发生变化时会对每个 $x_i$ 都产生影响，那么知道了 $df_1$ ，求解 $X$ 的一种解法是 $X$ 沿着 $X^{\prime}_{f_1}$ 的方向迭代就行了，公式（1）可以改写成：
$\begin{aligned} \to \left(X-X^{k}\right)=-X^{\prime}_{f_1}\sdot f\left(X^{k}\right) \\ \\ \to dX=-X^{\prime}_{f_1}\sdot f\left(X^{k}\right) \end{aligned}$
得到增量的表达式后，写成迭代的形式为：
$\begin{align} \to X^{k}=X^{k-1}+dX^{k-1} \end{align}$
其中
$\begin{align} dX^{k-1}=&-X^{\prime}_{f_1}\sdot f\left(X^{k-1}\right) \\ \nonumber\\ X^{\prime}_{f_1}=&P\odot\frac{1}{f^{\prime}\left(X^{k}\right)} \\ \nonumber\\ P=& \begin{bmatrix} 1/nx \\ 1/nx \\ \vdots \\ 1/nx \end{bmatrix}_{nx\times 1} \end{align}$
$\odot$ 是点乘， $P$ 的1-范数(每个元素绝对值之和)为1，其元素是对每个 $x_i$ 增量的权重，当 $X$ 沿着 $X^{\prime}_{f_1}$ 梯度的方向时均分，也可以根据迭代的效果设置。

2.2 求解偏导数矩阵

由于 $X$ 、 $f (X)$ 和 $\varphi(X)$ 均为列向量， $f (X)$ 对 $X$ 的偏导数写成：
$\begin{equation} f_{X}^{\prime}(X)= \begin{bmatrix} \frac{\partial f_{1}}{\partial x_{1}} & \frac{\partial f_{1}}{\partial x_{2}} & \cdots & \frac{\partial f_{1}}{\partial x_{nx}} \\ \frac{\partial f_{2}}{\partial x_{1}} & \frac{\partial f_{2}}{\partial x_{2}} & \cdots & \frac{\partial f_{ny}}{\partial x_{2}} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial f_{ny}}{\partial x_{1}} & \frac{\partial f_{ny}}{\partial x_{2}} & \cdots & \frac{\partial f_{ny}}{\partial x_{nx}} \end{bmatrix} \end{equation}$
使用 $f (x)$ 在 $x^{k}$ 的差分来近似函数在该点的导数，有几种差分：

后向差分
$\delta_{x}(x^{k})=\frac{f(x^{k})-f(x^{k}-\bigtriangleup x)}{\bigtriangleup x}$
前向差分
$\delta_{x}(x^{k})=\frac{f(x^{k}+\bigtriangleup x)-f(x^{k})}{\bigtriangleup x}$
中心差分
$\delta_{x}(x^{k})=\frac{f(x^{k}+\bigtriangleup x)-f(x^{k}-\bigtriangleup x)}{2\cdot \bigtriangleup x}$

2.3 步骤

步骤如下：

给定初值 $X^{0}=[x_{1}^{0},x_{2}^{0},\dots,x_{nx}^{0}]$ ，计算 $f(X^{0})=[f_{1}^{0},f_{2}^{0},\dots,f_{ny}^{0}]$ ；设置差分步长 $\bigtriangleup X=[\bigtriangleup x_{1},\bigtriangleup x_{2} ,\dots,\bigtriangleup x_{nx}]$ ，根据公式（8）计算偏导数 $f_{X}^{\prime}(X^{0})$ ；设置误差限 $\epsilon$ ，令初始迭代次数 $k = 1$ ;
根据公式（4）计算第 $k$ 次迭代 $X^{k}$ ，计算 $f(X^{k})$ ，根据公式（8）计算偏导数 $f_{X}^{\prime}(X^{k})$ ；
判断欧几里得范数 $\begin{Vmatrix} f(X^{k}) \end{Vmatrix}_{2}$ 是否小于 $\epsilon$ ，说人话就是 $f(X^{k})$ 是否等于0，如果满足，输出结果 $X^{k}$ ；如果不满足，则将迭代次数 $k$ 加1以更新迭代次数，重复步骤2和步骤3。
根据步骤1、步骤2所述的计算偏导数 $f_{X}^{\prime}(X^{k})$ ，可使用中心差分法来估计：
- 区别于 $X^{k}$ 表示第 $k$ 次迭代的初值，以下 $X^{nx+}$ 、 $X^{nx-}$ 表示对 $X^{k}$ 的第 $n x$ 个元素进行加、减增量， $f_{ny}^{nx+}$ 、 $f_{ny}^{nx-}$ 表示 $f(X^{nx+})$ 、 $f(X^{nx-})$ 的第 $n y$ 个元素。
- 增量
  
  $\begin{aligned} &X^{1+}=[x_{1}^{k}+\bigtriangleup x_{1},x_{2}^{k},\dots,x_{nx}^{k}]\to f(X^{1+})=[f_{1}^{1+},f_{2}^{1+},\dots,f_{ny}^{1+}] \\ \\ &X^{1-}=[x_{1}^{k}-\bigtriangleup x_{1},x_{2}^{k},\dots,x_{nx}^{k}]\to f(X^{1-})=[f_{1}^{1-},f_{2}^{1-},\dots,f_{ny}^{1-}] \\ \vdots \\ &X^{nx+}=[x_{1}^{k},x_{2}^{k},\dots,x_{nx}^{k}+\bigtriangleup x_{nx}]\to f(X^{nx+})=[f_{1}^{nx+},f_{2}^{nx+},\dots,f_{ny}^{nx+}] \\ \\ &X^{nx-}=[x_{1}^{k},x_{2}^{k},\dots,x_{nx}^{k}-\bigtriangleup x_{nx}]\to f(X^{nx-})=[f_{1}^{nx-},f_{2}^{nx-},\dots,f_{ny}^{nx-}] \end{aligned}$
- 中心差分
  
  $\begin{aligned} &\frac{\partial f_{j}}{\partial x_{1}} = \frac{f_{j}^{1+}-f_{j}^{1-}}{2\cdot \bigtriangleup x_{1}} \\ \\ &\frac{\partial f_{j}}{\partial x_{2}} = \frac{f_{j}^{2+}-f_{j}^{2-}}{2\cdot \bigtriangleup x_{2}} \\ &\vdots \\ &\frac{\partial f_{j}}{\partial x_{nx}}=\frac{f_{j}^{nx+}-f_{j}^{nx-}}{2\cdot \bigtriangleup x_{nx}} \\ \\ & j=1,\dots,ny \end{aligned}$
根据步骤1、步骤2所述的计算偏导数 $f_{X}^{\prime}(X^{k})$ ，亦可使用前向差分估计：
- 增量
  
  $\begin{aligned} &X^{1+}=[x_{1}^{k}+\bigtriangleup x_{1},x_{2}^{k},\dots,x_{nx}^{k}]\to f(X^{1+})=[f_{1}^{1+},f_{2}^{1+},\dots,f_{ny}^{1+}] \\ \\ &X^{2+}=[x_{1}^{k}+\bigtriangleup x_{1},x_{2}^{k}+\bigtriangleup x_{2},\dots,x_{nx}^{k}]\to f(X^{2+})=[f_{1}^{2+},f_{2}^{2+},\dots,f_{ny}^{2+}] \\ &\vdots \\ &X^{nx+}=[x_{1}^{k}+\bigtriangleup x_{1},x_{2}^{k}+\bigtriangleup x_{2},\dots,x_{nx}^{k}+\bigtriangleup x_{nx}]\to f(X^{nx+})=[f_{1}^{nx+},f_{2}^{nx+},\dots,f_{ny}^{nx+}] \end{aligned}$
- 前向差分
  
  $\begin{aligned} &\frac{\partial f_{j}}{\partial x_{1}} = \frac{f_{j}^{1+}-f_{j}^{k}}{ \bigtriangleup x_{1}} \\ \\ &\frac{\partial f_{j}}{\partial x_{2}} = \frac{f_{j}^{2+}-f_{j}^{1+}}{ \bigtriangleup x_{2}} \\ \vdots \\ &\frac{\partial f_{j}}{\partial x_{nx}} = \frac{f_{j}^{nx+}-f_{j}^{(nx-1)+}}{ \bigtriangleup x_{nx}} \\ \\ &j=1,\dots,ny \end{aligned}$

3. 实例测试

$\begin{aligned} f(x_1,x_2)&=(x_1-3)^{2}+(x_2-4)^{2}-25 \\ \\ X^0&=[-0.9,-0.9] \\ \\ \bigtriangleup X&=[0.0001, 0.0001] \\ \\ \epsilon&=0.001 \end{aligned}$
$P$ 取不同值时对应左右结果。


P=[0.5, 0.5]'	P=[1, 0]'

可见搜索方向不同得到不同满足要求的解。

4. 总结

介绍了牛顿-拉夫逊法在求解多输入单输出系统的实现。下回再讨论多输入多输出问题，以及一些应用场景。

5. 详细代码

ExampleNewtonLaphson.m

/*
* Auther:  sanfan66
*/
close all
clear all
%% 计算
x1=-1.2:0.1:1;
x2=-1.5:0.1:1;
[x,y]=meshgrid(x1,x2);
for i=1:size(x,1)
   for j=1:size(x,2) 
    f(i,j)=TestFunction([x(i,j),y(i,j)]');
   end
end

X0=[-0.9; -0.9];
deltaX0=[0.0001; 0.0001];
epsF=0.001;
try
    [index,X,fx,dx] = NewtonLaphsonMethod(X0,deltaX0,epsF,@TestFunction);
catch ME
     disp(ME.message);
end
%% 绘图
figure(1)
mesh(x,y,f)
hold on;
mesh(x,y,f*0)
hold on;
plot3(X(1,:),X(2,:),fx(:),'r-o')
hold on;
xCircle=-2:0.1:8;
plot3(xCircle,sqrt(25-(xCircle-3).^2)+4,xCircle*0,'r')
hold on;
plot3(xCircle,-sqrt(25-(xCircle-3).^2)+4,xCircle*0,'r')
xlabel("x1")
ylabel("x2")
zlabel("f")
xlim([-2 3])
ylim([-1 4])

figure
plot(1:index,fx(1:index),'-o')
xlabel("迭代次数")
ylabel("fx")
grid on;grid minor;
figure
plot(X(1,:),X(2,:),'-o')
xlabel("x1")
ylabel("x2")
% xlim([-0.46 -0.39])
% ylim([-0.34 -0.32])
grid on;grid minor;

fprintf("%5s%10s%10s%10s%10s\n","i","f0","ff","x1","x2");
for i=1:index
fprintf("%5.f%10f%10f%10f%10f\n",i,fx(1),fx(i),X(1,i),X(2,i));
end

function ft=TestFunction(x)
f=(x(1,:)-3).^2+(x(2,:)-4).^2-25;
ft=f';%转为列向量
end

NewtonLaphsonMethod.m

function [i,X,fx,dx] = NewtonLaphsonMethod(X0,deltaX0,epsF,objFun)
nx=size(X0,1);
xContribution=1/nx*ones(nx,1);
% xContribution=[1, 0]';
X(:,1)=X0;
fx(:,1)=objFun(X(:,1));
diffF(1,:)=SolveDiff(X0,deltaX0,objFun);
diffX(:,1)=1./diffF(1,:)';
dx(:,1)=diffX(:,1)*(-fx(:,1)).*xContribution;%将总变化量均分给两部分
maxIter=10;
for i=2:maxIter
    X(:,i)=X(:,i-1)+dx(:,i-1);
    fx(:,i)=objFun(X(:,i));
    diffF(i,:)=SolveDiff(X(:,i),deltaX0,objFun);
    diffX(:,i)=1./diffF(i,:)';
    dx(:,i)=diffX(:,i)*(-fx(:,i)).*xContribution;%将总变化量均分给两部分
    if(norm(fx(:,i))<=epsF)
        break;
    end
    
end

% if(i==maxIter && norm(fx(:,i))>epsF)
%     ME=MException('NewtonLaphsonError:maxIter', '超过最大迭代次数：i=%d，maxIter=%d',i+1,maxIter);
%     throw(ME);
% end

end

SolveDiff.m

function diffF = SolveDiff(X0,deltaX0,objFun)
%X0列向量
fx0=objFun(X0);
deltaX=diag(deltaX0);
XPlus=deltaX+repmat(X0,[1,length(X0)]);
XMinus=-deltaX+repmat(X0,[1,length(X0)]);
fxPlus=objFun(XPlus);
fxMinus=objFun(XMinus);
diffF=(fxPlus-fxMinus)./deltaX0/2;
end