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

文章目录

- 1. 概述
- 2. 原理步骤
- - 2.1 原理
  - 2.2 求解偏导数矩阵
  - 2.3 步骤
  - 3. 实例测试
  - 4. 总结
  - 5. 详细代码
    【非线性方程求解】割线法。
    1. 概述
    
    总结上篇博客内容，求解单变量非线性方程 f ( x ) = 0 f(x)=0 f(x)=0 时，牛顿法和割线法的区别在于当前导数的求解方法不同，之后都是用一条直线去近似 f ( x ) f(x) f(x) 。若求解的是非线性方程组，即含有多个输入和多个输出，实现有些不同。下面介绍牛顿法的具体实现，有时也叫作牛顿-拉夫逊（拉弗森）法。
    
    2. 原理步骤
    
    2.1 原理
    
    输入 { x i ∣ i = 1 , … , n x } \{x_{i}|i=1,\dots,nx\} {xi∣i=1,…,nx} 经过黑盒子后得到输出 { y j ∣ j = 1 , … , n y } \{y_{j}|j=1,\dots,ny\} {yj∣j=1,…,ny}
    
    但当我想要某种结果 { y j t ∣ j = 1 , … , n y } \{y_{j}^{t}|j=1,\dots,ny\} {yjt∣j=1,…,ny} 时，如何确定输入？为了解决这个问题，先按数学老师教的约定下符号：
    X = [ x 1 , … , x n x ] T X=[x_{1},\dots,x_{nx}]^{T} X=[x1,…,xnx]T
    T T T 为转秩符号，表示 X X X是一个列向量。同理：
    Y = [ y 1 , … , y n y ] T Y t = [ y 1 t , … , y n y t ] T \begin{aligned} Y=[y_{1},\dots,y_{ny}]^{T} \\ \\ Y^{t}=[y_{1}^{t},\dots,y_{ny}^{t}]^{T} \end{aligned} Y=[y1,…,yny]TYt=[y1t,…,ynyt]T
    t t t 表示目标输出。为了和之前讲过的步骤比较，把问题转换成求解以下方程，当然，你不想转换也行。
    f ( X ) = φ ( X ) − Y t = 0 f(X)=\varphi(X)-Y^{t}=0 f(X)=φ(X)−Yt=0
    接下来使用以前的老套路，仅仅把 x x x 换成 X X X 。为避免与 X X X 中的元素混淆，我把第 k k k 次迭代后的结果写成上标：
    → f ( X ) = f ( X k ) + f ′ ( X k ) ( X − X k ) + f ′ ′ ( ξ ) ( X − X k ) 2 / ( 2 ! ) → f ( X k ) + f ′ ( X k ) ( X − X k ) = 0 → f ′ ( X k ) ( X − X k ) = − f ( X k ) → X − X k = − f ( X k ) f ′ ( X k ) → X = X k − f ( X k ) f ′ ( 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} →f(X)=f(Xk)+f′(Xk)(X−Xk)+f′′(ξ)(X−Xk)2/(2!)→f(Xk)+f′(Xk)(X−Xk)=0→f′(Xk)(X−Xk)=−f(Xk)→X−Xk=−f′(Xk)f(Xk) →X=Xk−f′(Xk)f(Xk)
    
    你不能继续套路下去了，因为带有 X X X 的都是矩阵， f ′ ( X k ) f^{\prime}\left(X^{k}\right) f′(Xk)表示 f ( x ) f(x) f(x) 对 X X X 的偏导数矩阵，需要根据 f ′ ( X k ) f^{\prime}\left(X^{k}\right) f′(Xk) 的维度作不同变换，怎么检查 f ′ ( X k ) f^{\prime}\left(X^{k}\right) f′(Xk) 写的对不对呢？暂且将导数理解为分母变化量为1时分子的变化量，则当 X X X变化量为 d X dX dX 时 f ( X ) f(X) f(X)变化量为 d f ( X ) df(X) df(X)，检查下式是否成立即可。
    
    当系统为多输入单输出时，即 n x > 1 , n y = 1 nx>1,ny=1 nx>1,ny=1，如果这是线性方程，直接用线性代数的方法求解公式（1）就行了，仅需一步并且能证明有很多解。对于非线性方程，什么情况都有可能存在，需要迭代。公式（1）展开为：
    ∂ f 1 ∂ x 1 ⋅ d x 1 + ∂ f 1 ∂ x 2 ⋅ d x 2 + ⋯ + ∂ f 1 ∂ x n x ⋅ d x n x = d f 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} ∂x1∂f1⋅dx1+∂x2∂f1⋅dx2+⋯+∂xnx∂f1⋅dxnx=df1
    我们看到，每个 x i x_i xi 变化都会对 f 1 f_1 f1 产生影响，但问题我只想知道 f 1 = 0 f_1=0 f1=0 时 X = ？ X=？ X=？仿照公式（2）的形式，可以写出 X X X 对 f 1 f_1 f1 的导数矩阵 X f 1 ′ X^{\prime}_{f_1} Xf1′ 并且有：
    ∂ x 1 ∂ f 1 ⋅ d f 1 + ∂ x 2 ∂ f 1 ⋅ d f 1 + ⋯ + ∂ x n x ∂ f 1 ⋅ d f 1 = d X \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} ∂f1∂x1⋅df1+∂f1∂x2⋅df1+⋯+∂f1∂xnx⋅df1=dX
    公式(3)表示 f 1 f_1 f1 发生变化时会对每个 x i x_i xi 都产生影响，那么知道了 d f 1 df_1 df1 ，求解 X X X 的一种解法是 X X X 沿着 X f 1 ′ X^{\prime}_{f_1} Xf1′ 的方向迭代就行了，公式（1）可以改写成：
    → ( X − X k ) = − X f 1 ′ ⋅ f ( X k ) → d X = − X f 1 ′ ⋅ f ( X k ) \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} →(X−Xk)=−Xf1′⋅f(Xk)→dX=−Xf1′⋅f(Xk)
    得到增量的表达式后，写成迭代的形式为：
    → X k = X k − 1 + d X k − 1 \begin{align} \to X^{k}=X^{k-1}+dX^{k-1} \end{align} →Xk=Xk−1+dXk−1
    其中
    d X k − 1 = − X f 1 ′ ⋅ f ( X k − 1 ) X f 1 ′ = P ⊙ 1 f ′ ( X k ) P = [ 1 / n x 1 / n x ⋮ 1 / n x ] n x × 1 \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} dXk−1=Xf1′=P=−Xf1′⋅f(Xk−1)P⊙f′(Xk)1 1/nx1/nx⋮1/nx nx×1
    ⊙ \odot ⊙是点乘， P P P 的1-范数(每个元素绝对值之和)为1，其元素是对每个 x i x_i xi 增量的权重，当 X X X 沿着 X f 1 ′ X^{\prime}_{f_1} Xf1′ 梯度的方向时均分，也可以根据迭代的效果设置。
    
    2.2 求解偏导数矩阵
    
    由于 X X X 、 f ( X ) f(X) f(X) 和 φ ( X ) \varphi(X) φ(X)均为列向量， f ( X ) f(X) f(X) 对 X X X 的偏导数写成：
    f X ′ ( X ) = [ ∂ f 1 ∂ x 1 ∂ f 1 ∂ x 2 ⋯ ∂ f 1 ∂ x n x ∂ f 2 ∂ x 1 ∂ f 2 ∂ x 2 ⋯ ∂ f n y ∂ x 2 ⋮ ⋮ ⋱ ⋮ ∂ f n y ∂ x 1 ∂ f n y ∂ x 2 ⋯ ∂ f n y ∂ x n 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} fX′(X)= ∂x1∂f1∂x1∂f2⋮∂x1∂fny∂x2∂f1∂x2∂f2⋮∂x2∂fny⋯⋯⋱⋯∂xnx∂f1∂x2∂fny⋮∂xnx∂fny
    使用 f ( x ) f(x) f(x)在 x k x^{k} xk 的差分来近似函数在该点的导数，有几种差分：
    - 后向差分
      δ x ( x k ) = f ( x k ) − f ( x k − △ x ) △ x \delta_{x}(x^{k})=\frac{f(x^{k})-f(x^{k}-\bigtriangleup x)}{\bigtriangleup x} δx(xk)=△xf(xk)−f(xk−△x)
    - 前向差分
      δ x ( x k ) = f ( x k + △ x ) − f ( x k ) △ x \delta_{x}(x^{k})=\frac{f(x^{k}+\bigtriangleup x)-f(x^{k})}{\bigtriangleup x} δx(xk)=△xf(xk+△x)−f(xk)
    - 中心差分
      δ x ( x k ) = f ( x k + △ x ) − f ( x k − △ x ) 2 ⋅ △ x \delta_{x}(x^{k})=\frac{f(x^{k}+\bigtriangleup x)-f(x^{k}-\bigtriangleup x)}{2\cdot \bigtriangleup x} δx(xk)=2⋅△xf(xk+△x)−f(xk−△x)
      
      2.3 步骤
      
      步骤如下：
      1. 给定初值 X 0 = [ x 1 0 , x 2 0 , … , x n x 0 ] X^{0}=[x_{1}^{0},x_{2}^{0},\dots,x_{nx}^{0}] X0=[x10,x20,…,xnx0] ，计算 f ( X 0 ) = [ f 1 0 , f 2 0 , … , f n y 0 ] f(X^{0})=[f_{1}^{0},f_{2}^{0},\dots,f_{ny}^{0}] f(X0)=[f10,f20,…,fny0] ；设置差分步长 △ X = [ △ x 1 , △ x 2 , … , △ x n x ] \bigtriangleup X=[\bigtriangleup x_{1},\bigtriangleup x_{2} ,\dots,\bigtriangleup x_{nx}] △X=[△x1,△x2,…,△xnx] ，根据公式（8）计算偏导数 f X ′ ( X 0 ) f_{X}^{\prime}(X^{0}) fX′(X0) ；设置误差限 ϵ \epsilon ϵ ，令初始迭代次数 k = 1 k=1 k=1 ;
      2. 根据公式（4）计算第 k k k 次迭代 X k X^{k} Xk ，计算 f ( X k ) f(X^{k}) f(Xk) ，根据公式（8）计算偏导数 f X ′ ( X k ) f_{X}^{\prime}(X^{k}) fX′(Xk) ；
      3. 判断欧几里得范数 ∥ f ( X k ) ∥ 2 \begin{Vmatrix} f(X^{k}) \end{Vmatrix}_{2} f(Xk) 2 是否小于 ϵ \epsilon ϵ ，说人话就是 f ( X k ) f(X^{k}) f(Xk)是否等于0，如果满足，输出结果 X k X^{k} Xk ；如果不满足，则将迭代次数 k k k 加1以更新迭代次数，重复步骤2和步骤3。
      4. 根据步骤1、步骤2所述的计算偏导数 f X ′ ( X k ) f_{X}^{\prime}(X^{k}) fX′(Xk) ，可使用中心差分法来估计：
        
        区别于 X k X^{k} Xk表示第 k k k次迭代的初值，以下 X n x + X^{nx+} Xnx+、 X n x − X^{nx-} Xnx−表示对 X k X^{k} Xk的第 n x nx nx个元素进行加、减增量， f n y n x + f_{ny}^{nx+} fnynx+、 f n y n x − f_{ny}^{nx-} fnynx−表示 f ( X n x + ) f(X^{nx+}) f(Xnx+)、 f ( X n x − ) f(X^{nx-}) f(Xnx−)的第 n y ny ny个元素。
        
        增量
        
        X 1 + = [ x 1 k + △ x 1 , x 2 k , … , x n x k ] → f ( X 1 + ) = [ f 1 1 + , f 2 1 + , … , f n y 1 + ] X 1 − = [ x 1 k − △ x 1 , x 2 k , … , x n x k ] → f ( X 1 − ) = [ f 1 1 − , f 2 1 − , … , f n y 1 − ] ⋮ X n x + = [ x 1 k , x 2 k , … , x n x k + △ x n x ] → f ( X n x + ) = [ f 1 n x + , f 2 n x + , … , f n y n x + ] X n x − = [ x 1 k , x 2 k , … , x n x k − △ x n x ] → f ( X n x − ) = [ f 1 n x − , f 2 n x − , … , f n y n x − ] \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} ⋮X1+=[x1k+△x1,x2k,…,xnxk]→f(X1+)=[f11+,f21+,…,fny1+]X1−=[x1k−△x1,x2k,…,xnxk]→f(X1−)=[f11−,f21−,…,fny1−]Xnx+=[x1k,x2k,…,xnxk+△xnx]→f(Xnx+)=[f1nx+,f2nx+,…,fnynx+]Xnx−=[x1k,x2k,…,xnxk−△xnx]→f(Xnx−)=[f1nx−,f2nx−,…,fnynx−]
        
        中心差分
        
        ∂ f j ∂ x 1 = f j 1 + − f j 1 − 2 ⋅ △ x 1 ∂ f j ∂ x 2 = f j 2 + − f j 2 − 2 ⋅ △ x 2 ⋮ ∂ f j ∂ x n x = f j n x + − f j n x − 2 ⋅ △ x n x j = 1 , … , n y \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} ∂x1∂fj=2⋅△x1fj1+−fj1−∂x2∂fj=2⋅△x2fj2+−fj2−⋮∂xnx∂fj=2⋅△xnxfjnx+−fjnx−j=1,…,ny
        
        根据步骤1、步骤2所述的计算偏导数 f X ′ ( X k ) f_{X}^{\prime}(X^{k}) fX′(Xk) ，亦可使用前向差分估计：
        
        增量
        
        X 1 + = [ x 1 k + △ x 1 , x 2 k , … , x n x k ] → f ( X 1 + ) = [ f 1 1 + , f 2 1 + , … , f n y 1 + ] X 2 + = [ x 1 k + △ x 1 , x 2 k + △ x 2 , … , x n x k ] → f ( X 2 + ) = [ f 1 2 + , f 2 2 + , … , f n y 2 + ] ⋮ X n x + = [ x 1 k + △ x 1 , x 2 k + △ x 2 , … , x n x k + △ x n x ] → f ( X n x + ) = [ f 1 n x + , f 2 n x + , … , f n y n x + ] \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} X1+=[x1k+△x1,x2k,…,xnxk]→f(X1+)=[f11+,f21+,…,fny1+]X2+=[x1k+△x1,x2k+△x2,…,xnxk]→f(X2+)=[f12+,f22+,…,fny2+]⋮Xnx+=[x1k+△x1,x2k+△x2,…,xnxk+△xnx]→f(Xnx+)=[f1nx+,f2nx+,…,fnynx+]
        
        前向差分
        
        ∂ f j ∂ x 1 = f j 1 + − f j k △ x 1 ∂ f j ∂ x 2 = f j 2 + − f j 1 + △ x 2 ⋮ ∂ f j ∂ x n x = f j n x + − f j ( n x − 1 ) + △ x n x j = 1 , … , n y \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} ⋮∂x1∂fj=△x1fj1+−fjk∂x2∂fj=△x2fj2+−fj1+∂xnx∂fj=△xnxfjnx+−fj(nx−1)+j=1,…,ny
      3. 实例测试
      
      f ( x 1 , x 2 ) = ( x 1 − 3 ) 2 + ( x 2 − 4 ) 2 − 25 X 0 = [ − 0.9 , − 0.9 ] △ X = [ 0.0001 , 0.0001 ] ϵ = 0.001 \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} f(x1,x2)X0△Xϵ=(x1−3)2+(x2−4)2−25=[−0.9,−0.9]=[0.0001,0.0001]=0.001
      P P 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
```

导数 步骤

相关推荐

机器学习2---逻辑回归(基础准备）2024.03.10
PyQt5中下载Pyqt5-tools卡在Preparing metadata (pyproject.toml)-2024.03.09
微服务—ES数据同步2024.03.04
uniapp本地打包教程2024.03.03
如何在Python中安装OpenCV库2024.03.02

