【MATLAB学习笔记】数值方法——一维牛顿法（求极小值）

前言

  牛顿法是一种迭代求解极值问题的数值方法，它可以用来求解一维函数的极小值或极大值。该方法利用函数的一阶导数和二阶导数信息来逼近极值点。具体步骤如下：

1. 选择初始点 ${\alpha }_0$和允许误差$\epsilon$，并且令$k=0$；
2. 计算函数$f(x)$在${\alpha }_k$处的一阶导数$f^{\prime }({\alpha }_k)$和二阶导数$f^{\prime \prime }({\alpha }_k)$;
3. 求${\alpha }_{k+1}={\alpha }_k-\frac{f^{\prime }({\alpha }_k)}{f^{\prime \prime }({\alpha }_k)}$;
4. 若$\left|{\alpha }_{k+1}-{\alpha }_k\right|⩽\epsilon$，则可求得近似解${\alpha }^{\ast }={\alpha }_{k+1}$，停止计算，否则转至步骤5；
5. 令$k=k+1$，转回步骤2。

算法流程图

算法代码

function Newton_minValue(f_sym,x0,eps,N)
    % f_sym为所求函数, 符号函数
    % x0为初始点，默认为0
    % eps为允许误差, 默认为1e-4
    % N为最大迭代次数, 默认为500
    if nargin < 2, x0 = 0; end
    if nargin < 3, eps = 1e-4; end
    if nargin < 4, N = 500; end

    f = matlabFunction(f_sym);
    df = matlabFunction(diff(f_sym));
    ddf = matlabFunction(diff(f_sym,2));
    
    k = 0;
    fprintf('函数f=%s迭代开始:\n',f_sym)
    while true
        k = k + 1;

        x1 = x0 - df(x0)/(ddf(x0)+eps);
        fprintf('k=%2.d, x0=%.4f, df=%9.4f, ddf=%9.4f, x1=%.4f\n',k,x0,df(x0),ddf(x0),x1)
        if norm(abs(x1-x0),'fro')<eps
            fprintf('已收敛!结果如下：\nk=%2.d, x0=%.4f, x1=%.4f, f=%.4f',k,x0,x1,f(x1));
            break
        end

        if k > N
            fprintf('已经超出最大迭代次数%d！\n',N)
            break
        end
        x0 = x1;
    end
    fprintf('\n\n')

end

测试

测试函数如下：

syms x
f1_sym = x.^4 - 4*x.^3 - 6*x.^2 - 16*x + 4;
f2_sym = 4*x.^3 - 12*x.^2 - 12*x - 16;
f3_sym = x.^3 - 2*x + 1;
f4_sym = 3*x.^3 - 4*x + 5;

解析解

f1 = diff(f1_sym);
x1_min = double(vpa(solve(f1,x)));
f2 = diff(f2_sym);
x2_min = double(vpa(solve(f2,x)));
f3 = diff(f3_sym);
x3_min = double(vpa(solve(f3,x)));
f4 = diff(f4_sym);
x4_min = double(vpa(solve(f4,x)));

x1_min	x2_min	x3_min	x3_min
4.0000	2.4142	0.8165	0.6667

注：上表格仅展示了局部解。

算法解

x0 = 2;
eps = 1e-2;
N = 100;
Newton_minValue(f1_sym,x0,eps,N)
Newton_minValue(f2_sym,x0,eps,N)
Newton_minValue(f3_sym,x0,eps,N)
Newton_minValue(f4_sym,x0,eps,N)

迭代过程如下：

函数f=x^4 - 6*x^2 - 4*x^3 - 16*x + 4迭代开始:
k= 1, x0=2.0000, df= -56.0000, ddf= -12.0000, x1=-2.6706
k= 2, x0=-2.6706, df=-145.7204, ddf= 137.6760, x1=-1.6122
k= 3, x0=-1.6122, df= -44.6059, ddf=  57.8834, x1=-0.8417
k= 4, x0=-0.8417, df= -16.7867, ddf=  16.7034, x1=0.1627
k= 5, x0=0.1627, df= -18.2523, ddf= -15.5865, x1=-1.0091
k= 6, x0=-1.0091, df= -20.2209, ddf=  24.4388, x1=-0.1820
k= 7, x0=-0.1820, df= -14.2372, ddf=  -7.2331, x1=-2.1531
k= 8, x0=-2.1531, df= -85.7198, ddf=  95.3055, x1=-1.2538
k= 9, x0=-1.2538, df= -27.7020, ddf=  36.9546, x1=-0.5044
k=10, x0=-0.5044, df= -13.5134, ddf=   3.1574, x1=3.7621
k=11, x0=3.7621, df= -18.0024, ddf=  67.5479, x1=4.0285
k=12, x0=4.0285, df=   2.4266, ddf=  86.0645, x1=4.0003
k=13, x0=4.0003, df=   0.0291, ddf=  84.0249, x1=4.0000
已收敛!结果如下：
k=13, x0=4.0003, x1=4.0000, f=-156.0000

函数f=4*x^3 - 12*x^2 - 12*x - 16迭代开始:
k= 1, x0=2.0000, df= -12.0000, ddf=  24.0000, x1=2.4998
k= 2, x0=2.4998, df=   2.9925, ddf=  35.9950, x1=2.4167
k= 3, x0=2.4167, df=   0.0837, ddf=  34.0003, x1=2.4142
已收敛!结果如下：
k= 3, x0=2.4167, x1=2.4142, f=-58.6274

函数f=x^3 - 2*x + 1迭代开始:
k= 1, x0=2.0000, df=  10.0000, ddf=  12.0000, x1=1.1674
k= 2, x0=1.1674, df=   2.0882, ddf=   7.0042, x1=0.8696
k= 3, x0=0.8696, df=   0.2689, ddf=   5.2179, x1=0.8182
k= 4, x0=0.8182, df=   0.0084, ddf=   4.9093, x1=0.8165
已收敛!结果如下：
k= 4, x0=0.8182, x1=0.8165, f=-0.0887

函数f=3*x^3 - 4*x + 5迭代开始:
k= 1, x0=2.0000, df=  32.0000, ddf=  36.0000, x1=1.1114
k= 2, x0=1.1114, df=   7.1160, ddf=  20.0044, x1=0.7558
k= 3, x0=0.7558, df=   1.1413, ddf=  13.6046, x1=0.6720
k= 4, x0=0.6720, df=   0.0641, ddf=  12.0957, x1=0.6667
已收敛!结果如下：
k= 4, x0=0.6720, x1=0.6667, f=3.2222

总代码

clc;clear;close all
%% 
syms x
f1_sym = x.^4 - 4*x.^3 - 6*x.^2 - 16*x + 4;
f2_sym = 4*x.^3 - 12*x.^2 - 12*x - 16;
f3_sym = x.^3 - 2*x + 1;
f4_sym = 3*x.^3 - 4*x + 5;
% 解析解
f1 = diff(f1_sym);
x1_min = double(vpa(solve(f1,x)));
f2 = diff(f2_sym);
x2_min = double(vpa(solve(f2,x)));
f3 = diff(f3_sym);
x3_min = double(vpa(solve(f3,x)));
f4 = diff(f4_sym);
x4_min = double(vpa(solve(f4,x)));

x0 = 2;
eps = 1e-2;
N = 100;
Newton_minValue(f1_sym,x0,eps,N)
Newton_minValue(f2_sym,x0,eps,N)
Newton_minValue(f3_sym,x0,eps,N)
Newton_minValue(f4_sym,x0,eps,N)

%% 
function Newton_minValue(f_sym,x0,eps,N)
    % f_sym为所求函数, 符号函数
    % x0为初始点，默认为0
    % eps为允许误差, 默认为1e-4
    % N为最大迭代次数, 默认为500
    if nargin < 2, x0 = 0; end
    if nargin < 3, eps = 1e-4; end
    if nargin < 4, N = 500; end

    f = matlabFunction(f_sym);
    df = matlabFunction(diff(f_sym));
    ddf = matlabFunction(diff(f_sym,2));
    
    k = 0;
    fprintf('函数f=%s迭代开始:\n',f_sym)
    while true
        k = k + 1;

        x1 = x0 - df(x0)/(ddf(x0)+eps);
        fprintf('k=%2.d, x0=%.4f, df=%9.4f, ddf=%9.4f, x1=%.4f\n',k,x0,df(x0),ddf(x0),x1)
        if norm(abs(x1-x0),'fro')<eps
            fprintf('已收敛!结果如下：\nk=%2.d, x0=%.4f, x1=%.4f, f=%.4f',k,x0,x1,f(x1));
            break
        end

        if k > N
            fprintf('已经超出最大迭代次数%d！\n',N)
            break
        end
        x0 = x1;
    end
    fprintf('\n\n')

end

总结

  这只是一个基础的示例，实际中还会有更具体的、更细致的要求，这就需要再做额外调整；另外本人也仍在学习中，这只是个人的学习笔记，可能还有一些不足之处，欢迎指正。

参考文献

主编白清顺. 机械优化设计.第6版 [M]. 2017.