目录
1/3
2/3
3/3
MATLAB程序设计——牛顿法求解非线性方程
鸣凤在竹-白驹食场
2022年04月30日 16:25
收录于文集
共6篇

1 牛顿迭代法

牛顿迭代法实质上是一种线性化方法，其基本思想是将非线性方程f(x)%3D0逐步归结为某种线性方程来求解。

1.1 牛顿法

牛顿迭代法又称切线法，是一种有特色的求根方法。用牛顿迭代法求f(x)%3D0的单根x%5E%7B*%7D的主要步骤：

（1）Newton法的迭代公式

x_%7Bk%2B1%7D%3D%20x_%7Bk%7D-%5Cfrac%7Bf(x_%7Bk%7D)%7D%7Bf%26%2339%3B(x_%7Bk%7D)%20%7D%20%EF%BC%8Ck%3D0%2C1%2C2%2C...

（2）以x%5E%7B*%7D附近的某一个值x_0为迭代初值，代入迭代公式，反复迭代，得到序列x_1%2C%20x_2%2C%20x_3%2C...

（3）若序列收敛，则必收敛于精确根x%5E%7B*%7D，即%5Clim_%7Bk%5Cto%5Cinfty%20%7D%20x_k%20%3Dx%5E*。

牛顿法有显然的几何意义，方程f(x)%3D0的根x%5E%7B*%7D可解释为曲线与轴交点的横坐标。设x_k是根x%5E%7B*%7D的某个近似值，过曲线y%3Df(x)上横坐标为x_k的点P_k(x_k%2C%20y_k)引切线，并将该切线与x轴交点的横坐标x_%7Bk%2B1%7D作为x%5E%7B*%7D的新的近似值。

牛顿法初值的选择：若f(x)在%5Ba%2Cb%5D上连续，存在2阶导数，且满足下列条件：

（1）f(a)f(b)%3C0；

（2）f%26%2339%3B(x)在(a%2Cb)内不变号，且f%26%2339%3B(x)%5Cneq%200；

（3）f%26%2339%3B%26%2339%3B(x)在(a%2Cb)内不变号，且f%26%2339%3B%26%2339%3B(x)%5Cneq%200；

（4）%5Cfrac%7B%5Cvert%20f(a)%20%5Cvert%20%7D%7B%5Cvert%20f%26%2339%3B(a)%20%5Cvert%20%7D%20%5Cleq%20b-a，且 %5Cfrac%7B%5Cvert%20f(b)%20%5Cvert%20%7D%7B%5Cvert%20f%26%2339%3B(b)%20%5Cvert%20%7D%20%5Cleq%20b-a；

则对任意的初值x_0%5Cin%20%5Ba%2Cb%5D，牛顿迭代序列收敛于f(x)%3D0在%5Ba%2Cb%5D内的唯一根。

牛顿的优缺点：牛顿法是目前求解非线性方程（组）的主要方法，至少二阶局部收敛，收敛速度较快，特别是当迭代初值充分靠近精确解时。对重根收敛速度较慢（线性收敛），对初值的选取很敏感，要求初值相当接近真解，需要求导数。

图1 牛顿法迭代收敛过程的几何意义

1.2 哈利牛顿法

哈利法（Halley）可用于加速牛顿法收敛，哈利迭代公式：

x_%7Bk%2B1%7D%3D%20x_%7Bk%7D-%5Cfrac%7Bf(x_%7Bk%7D)%7D%7Bf%26%2339%3B(x_%7Bk%7D)%20%7D%20%5Cleft(1-%5Cfrac%7Bf(x_k)f%26%2339%3B%26%2339%3B(x_k)%7D%7B2(f%26%2339%3B(x_k))%5E2%7D%20%5Cright)%5E%7B-1%7D%EF%BC%8Ck%3D0%2C1%2C2%2C...

哈利法在单根情况下可达到三阶收敛。

1.3 简化牛顿法

将牛顿法的迭代公式改为

x_%7Bk%2B1%7D%3Dx_k-%5Clambda%20f(x_k)%EF%BC%8Ck%3D0%2C1%2C2%2C...

为保证收敛，系数%5Clambda%20只需满足0%3C%5Clambda%20%3C%5Cfrac%7B2%7D%7Bf%26%2339%3B(x_k)%7D%20即可。如果%5Clambda%20取常数%5Cfrac%7B1%7D%7Bf%26%2339%3B(x_0)%7D%20，则称为简化牛顿法，也称平行弦法。

简化牛顿法线性收敛。

1.4 牛顿下山法

牛顿法的收敛性依赖初值x_0的选取，如果x_0偏离所求根x%5E*较远，则牛顿法可能发散。牛顿下山法，引入下山因子%5Clambda%20，保证函数值稳定下降的同时，加速收敛速度。

牛顿下山法公式：

x_%7Bk%2B1%7D%3D%20x_%7Bk%7D-%5Clambda%20%5Cfrac%7Bf(x_%7Bk%7D)%7D%7Bf%26%2339%3B(x_%7Bk%7D)%20%7D%20%EF%BC%8Ck%3D0%2C1%2C2%2C...

下山因子的取法：从%5Clambda%20%3D1开始，逐次减半，即%5Clambda%20%3D1%2C%5Cfrac%7B1%7D%7B2%7D%20%2C%5Cfrac%7B1%7D%7B2%5E2%7D%20%2C%5Cfrac%7B1%7D%7B2%5E3%7D%20%2C...，直到满足下降条件。

1.5 重根情形

当x%5E*为f(x)%3D0的m(m%3E0)重根时，则f(x)可表为f(x)%3D(x-x%5E*)%5Emg(x)，其中g(x%5E*)%5Cneq%200，此时用牛顿迭代法求x%5E*仍然收敛，只是收敛速度将大大减慢，牛顿法求方程的重根时仅为线性收敛。

将求重根问题化为求单根问题进行牛顿法求解。对f(x)%3D(x-x%5E*)%5Emg(x)，令函数

%5Cmu%20(x)%3D%5Cfrac%7Bf(x)%7D%7Bf%26%2339%3B(x)%7D%3D%5Cfrac%7B(x-x%5E*)g(x)%7D%7Bmg(x)%2B(x-x%5E*)g%26%2339%3B(x)%7D%20%20

则化为求%5Cmu%20(x)%3D0的单根x%5E*的问题，对它用牛顿法是二阶收敛的。其迭代函数为

%5Cvarphi%20(x)%3Dx-%5Cfrac%7B%5Cmu(x)%20%7D%7B%5Cmu%26%2339%3B(x)%7D%20%3Dx-%5Cfrac%7Bf(x)f%26%2339%3B(x)%7D%7B%5Bf%26%2339%3B(x)%5D%5E2-f(x)f%26%2339%3B%26%2339%3B(x)%7D%20

从而构造出迭代方法

x_%7Bk%2B1%7D%3Dx_k-%5Cfrac%7Bf(x_k)f%26%2339%3B(x_k)%7D%7B%5Bf%26%2339%3B(x_k)%5D%5E2-f(x_k)f%26%2339%3B%26%2339%3B(x_k)%7D%EF%BC%8C%20k%3D0%2C1%2C2%2C...

2 牛顿迭代法MATLAB算法实现

算法说明：

（1）输入参数的判定和必要的处理，包括可变参数的处理，若某个可变参数不输入，则采用默认值。

（2）由于牛顿法需要计算方程的一阶导数和二阶导数（重根情形），故方程的定义采用符号定义（避免手工计算方程的导数所带来的计算量），在算法solve_diff_fun(equ)中实现符号函数的一阶导数和二阶导数的计算，并转换为匿名函数进行数值运算。

（3）根据参数method选择执行不同的牛顿迭代方法，具体：

    1）“newton”为牛顿法newton()；

    2）“simplify”为简化的牛顿法simplify_newton()；

    3）“halley”为哈利法（牛顿加速法）newton_halley()；

    4）“downhill”为牛顿下山法newton_downhill()，下山法包括了下山因子的存储；

    5）“multi”为重根情形newton_multi_root()。

（4）根据参数display选择是否显示迭代过程信息或仅显示结果信息。

（5）可变输出参数的判定和必要处理。

思考：MATLAB中没有类似于Python的list结构，如何在未知迭代次数的情况下，预分配内存，存储迭代过程变量信息？

代码块
JavaScript
自动换行
复制代码
function  varargout = newton_root(equ_func, x0, varargin)
%% 牛顿法求解方程的根，包含简化牛顿法 simplify ，牛顿法 newton ，牛顿加速哈利法 halley ，牛顿下山法 downhill 和重根情形 multi_r 。
% 1. equ_func表示待求（非）线性方程，要求是符号函数定义
% 2. x0表示迭代求解的初值
% 3. varargin  表示可变参数
% 4. 输出参数varargout为可变参数：
% 4.1 如果输出参数为3个，则第1个参数为近似解，第2个参数为近似解的精度误差，第3个参数为退出条件
% 4.2 如果输出参数为2个，则第1个参数为近似解，第2个参数为近似解的精度误差
% 4.3 如果输出参数为1个，则表示牛顿法求解的结构体，包括各种迭代过程中的信息

%% 1.输入参数的判别与处理
if nargin < 2
    error("输入参数不足，请输入完整参数equ_func和x0....")
end

if class(equ_func) ~= "sym"
    error("缺失参数 equ_func 或参数 equ_func 定义错误...")
end

if class(x0) ~= "double"
    error("缺失参数 x0 或参数 x0 定义错误...")
end

args = inputParser; % 函数的输入解析器
addParameter(args, 'method', "newton"); % 设置牛顿法的求解形式，默认值普通newton法
addParameter(args, 'eps', 1e-10); % 设置求解精度参数eps，默认值1e-10
addParameter(args, 'display', 'final'); % 设置是否显示迭代过程参数display，默认值final，可为iter
addParameter(args, 'max_iter', 200); % 设置最大迭代次数参数max_iter，默认值200
parse(args,varargin{:}); % 对输入变量进行解析，如果检测到前面的变量被赋值，则更新变量取值

method = args.Results.method;  % 求解方法
global  eps  % 精度，定义为全局变量
eps = args.Results.eps;  % 若未传参，精度获取默认值
display = args.Results.display;  % 是否显示迭代过程
global max_iter  % 最大迭代次数，全局变量
max_iter = int16(args.Results.max_iter);  % 若未传参，最大迭代次数获取默认值

%% 2. 计算f(x)的一阶导和二阶导函数，并转换为匿名函数
[fx, dfx, d2fx] = solve_diff_fun(equ_func);

%% 3. 根据所采用的求解方法method选择对应的牛顿法形式求解
if lower(method) == "newton"
    method_message = "普通牛顿迭代法";
    iter_info = newton(fx, dfx, x0);  % 普通牛顿法
elseif lower(method) == "halley"
    method_message = "牛顿—哈利迭代法";
    iter_info = newton_halley(fx, dfx, d2fx, x0);  % 牛顿—哈利法
elseif lower(method) == "simplify"
    method_message = "简化的牛顿迭代法";
    iter_info = simplify_newton(fx, dfx, x0);  % 简化牛顿法
elseif lower(method) == "downhill"
    method_message = "牛顿下山迭代法";
    [iter_info, downhill_lambda] = newton_downhill(fx, dfx, x0);  % 牛顿下山法
elseif lower(method) == "multi_r"
    method_message = "牛顿重根情形迭代法";
    iter_info = newton_multi_root(fx, dfx, d2fx, x0);  % 牛顿重根情形
else
    warning("牛顿法仅限于【newton、halley、simple、downhill和multi_r】之一，此处采用牛顿法求解！")
    method_message = "普通牛顿迭代法";
    iter_info = newton(fx, dfx, x0);  % 普通牛顿法
end

%% 4. 根据参数display，显示迭代过程或迭代结果信息
if lower(display) == "iter"
    disp(method_message)
    disp(array2table(iter_info, 'VariableNames', ["IterNums", "ApproximateSol", "ErrorPrecision"]))
elseif lower(display) == "final"
    fprintf(method_message + "：IterNum = %d,  ApproximateSolution = %.15f,  ErrorPrecision = %.15e\n", ...
        iter_info(end, 1), iter_info(end, 2), iter_info(end, 3));
else
    warning("参数display仅限于【iter、final】之一！")
end

%% 2. 1求符号方程的一阶导和二阶导，并转换为匿名函数，方便计算
    function [equ_expr, diff_equ, diff2_equ] = solve_diff_fun(equ)
        x = symvar(equ);    % 获取符号函数的变量
        diff_equ = matlabFunction(diff(equ, x, 1));  % 一阶导数方程，并转换为匿名函数
        diff2_equ = matlabFunction(diff(equ, x, 2));  % 二阶导数方程，并转换为匿名函数
        equ_expr = matlabFunction(equ);  % 原符号方程转换为匿名函数
    end

%% 3. 1 牛顿法求解
    function iter_info = newton(fx, dfx, x0)
        var_iter = 0;  % 迭代次数变量
        sol_tol = fx(x0);  % 方程在当前迭代值的精度，初始为初值时的精度
        if abs(sol_tol) < eps
            iter_info = [1, x0, sol_tol];
            return
        end
        x_cur = x0;  % 迭代过程中的对应于x_(k)
        iter_info =zeros(max_iter, 3);  % 存储迭代过程信息，思考如何预分配内存？
        
        % 算法具体迭代多少次，不清楚，故用while。可用for循环，在循环体内，满足精度要求退出循环
        while abs(sol_tol) > eps && var_iter < max_iter
            var_iter = var_iter + 1;  % 迭代次数增一
            x_next = x_cur - fx(x_cur) / dfx(x_cur);  % 牛顿迭代公式
            sol_tol = fx(x_next);  % 新值的精度
            x_cur = x_next;  % 近似根的迭代更替
            iter_info(var_iter, :) = [var_iter, x_next, sol_tol];
        end
        
        if var_iter < max_iter
            iter_info(var_iter + 1:end, :) = [];  % 删除未赋值的区域
        end
    end

%% 3. 2 牛顿—哈利法法求解
    function iter_info = newton_halley(fx, dfx, d2fx, x0)
        var_iter = 0;  % 迭代次数变量
        sol_tol = fx(x0);  % 方程在当前迭代值的精度，初始为初值时的精度
        if abs(sol_tol) < eps
            iter_info = [1, x0, sol_tol];
            return
        end
        x_cur = x0;  % 迭代过程中的对应于x_(k)
        iter_info =zeros(max_iter, 3);  % 存储迭代过程信息
        
        while abs(sol_tol) > eps && var_iter < max_iter
            var_iter = var_iter + 1;  % 迭代次数增一
            fval = fx(x_cur);  % 当前近似解的方程精度
            df_val = dfx(x_cur);  % 当前近似解的方程一阶导数精度
            d2f_val = d2fx(x_cur);  % 当前近似解的方程二阶导数精度
            x_next = x_cur - fval / df_val / (1 - fval * d2f_val / (2* df_val^2));  % 牛顿哈利迭代公式
            sol_tol = fx(x_next);  % 新值的精度
            x_cur = x_next;  % 近似根的迭代更替
            iter_info(var_iter, :) = [var_iter, x_next, sol_tol];
        end
        
        if var_iter < max_iter
            iter_info(var_iter + 1:end, :) = [];  % 删除未赋值的区域
        end
    end

%% 3. 3 简化牛顿法求解
    function iter_info = simplify_newton(fx, dfx, x0)
        var_iter = 0;  % 迭代次数变量
        lambda = 1 / dfx(x0);  % 简化牛顿法的常数
        sol_tol = fx(x0);  % 方程在当前迭代值的精度，初始为初值时的精度
        if abs(sol_tol) < eps
            iter_info = [1, x0, sol_tol];
            return
        end
        x_cur = x0;  % 迭代过程中的对应于x_(k)
        iter_info =zeros(max_iter, 3);  % 存储迭代过程信息，思考如何预分配内存？
        
        while abs(sol_tol) > eps && var_iter < max_iter
            var_iter = var_iter + 1;  % 迭代次数增一
            x_next = x_cur - lambda * fx(x_cur);  % 简化牛顿迭代公式
            sol_tol = fx(x_next);  % 新值的精度
            x_cur = x_next;  % 近似根的迭代更替
            iter_info(var_iter, :) = [var_iter, x_next, sol_tol];
        end
        
        if var_iter < max_iter
            iter_info(var_iter + 1:end, :) = [];  % 删除未赋值的区域
        end
    end

%% 3. 4 牛顿下山法求解
    function [iter_info, downhill_lambda] = newton_downhill(fx, dfx, x0)
        var_iter = 0;  % 迭代次数变量
        sol_tol = fx(x0);  % 方程在当前迭代值的精度，初始为初值时的精度
        if abs(sol_tol) < eps
            iter_info = [1, x0, sol_tol];
            downhill_lambda = [1];
            return
        end
        x_cur = x0;  % 迭代过程中的对应于x_(k)
        iter_info =zeros(max_iter, 3);  % 存储迭代过程信息，思考如何预分配内存？
        downhill_lambda = ones(max_iter);  % 存储下山因子
        
        while abs(sol_tol) > eps && var_iter < max_iter
            var_iter = var_iter + 1;  % 迭代次数增一
            f_val = fx(x_cur);  % 当前近似解的方程精度
            df_val = dfx(x_cur);  % 当前近似解的方程一阶导数精度
            x_next = x_cur - f_val / df_val;  % 牛顿迭代公式
            sol_tol = fx(x_next);  % 新值的精度
            
            % 判别是否满足下降条件
            while abs(sol_tol) > abs(f_val)
                downhill_lambda(var_iter) = downhill_lambda(var_iter) / 2;  % 下山因子减半
                x_next = x_cur - downhill_lambda(var_iter) * f_val / df_val;  % 牛顿下山迭代公式
                sol_tol = fx(x_next);  % 新值的精度
            end
            
            x_cur = x_next;  % 近似根的迭代更替
            iter_info(var_iter, :) = [var_iter, x_next, sol_tol];
        end
        
        if var_iter < max_iter
            iter_info(var_iter + 1:end, :) = [];  % 删除未赋值的区域
            downhill_lambda(var_iter + 1:end) = [];  % 删除未赋值的区域
        end
    end

%% 3. 5 牛顿法重根情形求解
    function iter_info = newton_multi_root(fx, dfx,d2fx, x0)
        var_iter = 0;  % 迭代次数变量
        sol_tol = fx(x0);  % 方程在当前迭代值的精度，初始为初值时的精度
        if abs(sol_tol) < eps
            iter_info = [1, x0, sol_tol];
            return
        end
        x_cur = x0;  % 迭代过程中的对应于x_(k)
        iter_info =zeros(max_iter, 3);  % 存储迭代过程信息，思考如何预分配内存？
        
        while abs(sol_tol) > eps && var_iter < max_iter
            var_iter = var_iter + 1;  % 迭代次数增一
            fval = fx(x_cur);  % 当前近似解的方程精度
            df_val = dfx(x_cur);  % 当前近似解的方程一阶导数精度
            d2f_val = d2fx(x_cur);  % 当前近似解的方程二阶导数精度
            x_next = x_cur -fval * df_val / (df_val^2 - fval * d2f_val);  % 牛顿重根情形迭代公式
            sol_tol = fx(x_next);  % 新值的精度
            x_cur = x_next;  % 近似根的迭代更替
            iter_info(var_iter, :) = [var_iter, x_next, sol_tol];
        end
        
        if var_iter < max_iter
            iter_info(var_iter + 1:end, :) = [];  % 删除未赋值的区域
        end
    end

%% 5. 输出参数结构体的构造
sol.ApproximateSolution = iter_info(end, 2);  % 数值近似解
sol.ErrorPrecision =iter_info(end, 3);  % 数值近似解误差
sol.NumerOfIterations = iter_info(end, 1);  % 迭代次数
sol.NewtonMethod = method;  % 所采用的方法
sol.IntrationProcessInfomation = iter_info;  % 迭代过程信息
if length(iter_info) == max_iter
    sol.info = "已达最大迭代次数";
end
if length(iter_info) < max_iter
    sol.convergence = "收敛到满足精度的近似解";
else
    sol.convergence = "精度过高，或未收敛到满足精度的解，增大迭代次数或修改初值或降低精度要求";
end
if lower(method) == "downhill"
    sol.DownhillLambda = downhill_lambda;
end

%% 6. 输出参数的判别
if nargout == 3
    varargout{1} = iter_info(end, 2);
    varargout{2} = iter_info(end, 3);
    if length(iter_info) == max_iter
        varargout{3} = 0;  % 退出标记
    else
        varargout{3} = 1; % 退出标记
    end
elseif nargout == 2
    varargout{1} = iter_info(end, 2);
    varargout{2} = iter_info(end, 3);
elseif nargout == 1
    varargout{1} = sol;
elseif nargout == 0
	 varargout{1} = iter_info(end, 2); % 仅给出近似解
else
    error("输出参数最多为3个....")
end
end
复制成功

3 牛顿迭代法案例测试

例1：求非线性方程在x%5Cin%20%5B-3.5%2C%205%5D区间的近似根

f(x)%3D2e%5E%7B-x%7Dsinx%2B2cosx-0.25%3D0

可视化函数图像：

代码块
JavaScript
自动换行
复制代码
syms x
equ_func = 2 * exp(-x) * sin(x) + 2 * cos(x) - 0.25;
fplot(equ_func, [-3.5, 5], "r-", "LineWidth", 1.5);
hold on
fplot(@(x) x * 0, [-3.5, 5], "k--")
xlabel('$x$','FontSize',12,'Interpreter','latex');
ylabel('$y$','FontSize',12,'Interpreter','latex'); 
title("$ 2e^{-x}sinx + 2cosx - 0.25=0$", 'FontSize', 14, 'Interpreter','latex')
grid on; grid minor; box off; hold off
复制成功

图2 非线性方程图像

从图中可以看出，在区间内存在四个根，此处仅求在x_0%3D0附近的根，其他根求解自行设计。

验证输入参数：

代码块
JavaScript
自动换行
复制代码
>> sol_newton = newton_root(equ_func, "display", "iter", "eps", 1e-16, "method", "newton");
错误使用 newton_root (第 21 行)
缺失参数 x0 或参数 x0 定义错误... 
>> sol_newton = newton_root(x0, "display", "iter", "eps", 1e-16, "method", "newton");
错误使用 newton_root (第 17 行)
缺失参数 equ_func 或参数 equ_func 定义错误...
复制成功

验证输出参数，仅有一个输出参数，则为结构体，且输出迭代过程：

代码块
JavaScript
自动换行
复制代码
>> sol_newton = newton_root(equ_func, x0, "display", "iter", "eps", 1e-16, "method", "newton")
普通牛顿迭代法
    IterNums      ApproximateSol         ErrorPrecision    
    ________    __________________    _____________________
       1                    -0.875        -2.65048857279132
       2        -0.555391102428116       -0.388296167451843
       3        -0.489056835053933      -0.0167045381738624
       4        -0.485935114436592    -3.66522867834007e-05
       5        -0.485928234722167    -1.77926784417082e-10
       6         -0.48592823468877                        0
sol_newton = 
  包含以下字段的 struct:

           ApproximateSolution: -0.485928234688770
                ErrorPrecision: 0
             NumerOfIterations: 6
                  NewtonMethod: "newton"
    IntrationProcessInfomation: [6×3 double]
                   convergence: "收敛到满足精度的近似解"
复制成功

验证输出参数，有三个输出参数，分别为近似解、方程在近似解处的误差精度和退出标记，1表示收敛到满足精度要求的近似解，且输出最终迭代结果：

代码块
JavaScript
自动换行
复制代码
>> [x, f_eps, exitflag, iter] = newton_root(equ_func, x0, "display", "final", "eps", 1e-16, "method", "newton")
普通牛顿迭代法：IterNum = 6,  ApproximateSolution = -0.485928234688770,  ErrorPrecision = 0.000000000000000e+00
错误使用 newton_root (第 235 行)
输出参数最多为3个.... 
>> [x, f_eps, exitflag] = newton_root(equ_func, x0, "display", "final", "eps", 1e-16, "method", "newton")
普通牛顿迭代法：IterNum = 6,  ApproximateSolution = -0.485928234688770,  ErrorPrecision = 0.000000000000000e+00
x =
  -0.485928234688770
f_eps =
     0
exitflag =
     1
复制成功

例2：求带有重根的非线性方程在x_0%3D0.9附近的近似解

f(x)%3D(x-1)*(sin(x-1)%2B3x)-x%5E3%2B1%3D0

由于简化的牛顿法收敛速度较慢，此处不再可视化简化的牛顿法。

代码块
JavaScript
自动换行
复制代码
syms x
equ_func = (x - 1) * (sin(x - 1) + 3 * x) - x^3 + 1;
x0 = 0.9;
sol_newton = newton_root(equ_func, x0, "display", "iter", "eps", 1e-16, "method", "newton");
sol_halley = newton_root(equ_func, x0, "display", "iter", "eps", 1e-16, "method", "halley");
sol_simplify = newton_root(equ_func, x0, "display", "iter", "eps", 1e-16, "method", "simplify");
sol_downhill = newton_root(equ_func, x0, "display", "iter", "eps", 1e-16, "method", "downhill");
sol_multi_r = newton_root(equ_func, x0, "display", "iter", "eps", 1e-16, "method", "multi_r");
figure
plot(sol_newton.IntrationProcessInfomation(:, 1), sol_newton.IntrationProcessInfomation(:, 3), "o-", "LineWidth", 1)
hold on
plot(sol_halley.IntrationProcessInfomation(:, 1), sol_halley.IntrationProcessInfomation(:, 3), "s-", "LineWidth", 1)
% plot(sol_simplify.IntrationProcessInfomation(:, 1), sol_simplify.IntrationProcessInfomation(:, 3), "p-", "LineWidth", 1)
plot(sol_downhill.IntrationProcessInfomation(:, 1), sol_downhill.IntrationProcessInfomation(:, 3), "p-", "LineWidth", 1)
plot(sol_multi_r.IntrationProcessInfomation(:, 1), sol_multi_r.IntrationProcessInfomation(:, 3), "*-", "LineWidth", 1)
legend("newton：" + string(sol_newton.NumerOfIterations), "halley：" + string(sol_halley.NumerOfIterations), ...
    "downhill：" + string(sol_downhill.NumerOfIterations), "multiroot：" + string(sol_multi_r.NumerOfIterations))
legend("boxoff")
grid on; grid minor
box off; hold off
axis([0, 25, -0.0003, 3e-3])
xlabel("Numer Of Iterations")
ylabel("Error Precision")
title("The Different Newton Method of Error Precision Curve")
复制成功

可视化各牛顿法的误差精度下降曲线如图3，可见在存在重根情形下，采用牛顿重根公式，收敛速度极快，其次为牛顿哈利法。由于此处未用到下山因子，故牛顿法和牛顿下山法求解过程一致。

图3 带有重根的牛顿法各迭代方法的精度误差曲线

各方法迭代过程

代码块
JavaScript
自动换行
复制代码
普通牛顿迭代法
    IterNums     ApproximateSol         ErrorPrecision   
    ________    _________________    ____________________
        1       0.947892373805653     0.00285545902358852
        2       0.973327008365232    0.000730340585354261
        3       0.986493966243405    0.000184871077568705
        4       0.993202483991711    4.65199555944595e-05
        5       0.996589832769599    1.16688756544026e-05
        6       0.998292027147433    2.92215229613646e-06
        7       0.999145286557988    7.31159377154178e-07
        8          0.999572460931    1.82867799680686e-07
        9       0.999786184803889    4.57267127496053e-08
       10       0.999893080977268    1.14328996270174e-08
       11       0.999946537631102    2.85837764568697e-09
       12        0.99997326810061    7.14613479502191e-10
       13       0.999986633870396    1.78655867877353e-10
       14       0.999993316892757    4.46641612583676e-11
       15        0.99999665843075      1.116617909247e-11
       16       0.999998329221071    2.79143375081503e-12
       17       0.999999164587949    6.97886193279373e-13
       18       0.999999582277231    1.74527059471075e-13
       19       0.999999791180077    4.36317648677687e-14
       20       0.999999895652271    1.08801856413265e-14
       21       0.999999947786534    2.66453525910038e-15
       22       0.999999973302319     7.7715611723761e-16
       23       0.999999987857067    1.11022302462516e-16
       24       0.999999992428546                       0
牛顿—哈利迭代法
    IterNums     ApproximateSol         ErrorPrecision   
    ________    _________________    ____________________
        1       0.965551810226976     0.00122732199318976
        2       0.988384456906306    0.000136484986290042
        3       0.996113110670598    1.51665933865175e-05
        4       0.998702689489662    1.68519748056095e-06
        5       0.999567376081915    1.87244420080113e-07
        6       0.999855771228448    2.08049387717679e-08
        7       0.999951921431764    2.31165986352977e-09
        8       0.999983973553835    2.56851095947752e-10
        9       0.999994657822902    2.85390600041069e-11
       10       0.999998219279682    3.17101900293437e-12
       11       0.999999406450421    3.52273765713562e-13
       12        0.99999980210885     3.9190872769268e-14
       13       0.999999934170877    4.32986979603811e-15
       14       0.999999978008308    4.44089209850063e-16
       15       0.999999991113451    1.11022302462516e-16
       16        1.00000000074541                       0
简化的牛顿迭代法
    IterNums     ApproximateSol         ErrorPrecision   
    ________    _________________    ____________________
        1       0.947892373805653     0.00285545902358852
        2       0.960343476131293     0.00163459314954228
        3       0.967471047004351     0.00109236613692898
        4       0.972234261057016    0.000792242820659372
        5       0.975688800710815    0.000605344947669506
        6       0.978328380485372    0.000479800603124803
        7       0.980420529718144    0.000390837064078497
        8       0.982124757387774    0.000325218857110121
        9       0.983542859778547    0.000275282449072156
       10       0.984743216678623    0.000236311720020921
       ......
      196       0.998877331341468    1.26179964665685e-06
      197       0.998882833362819    1.24945532808951e-06
      198       0.998888281557317     1.2372916339265e-06
      199       0.998893676712566    1.22530505242135e-06
      200       0.998899019600859    1.21349215642663e-06
牛顿下山迭代法
    IterNums     ApproximateSol         ErrorPrecision   
    ________    _________________    ____________________
        1       0.947892373805653     0.00285545902358852
        2       0.973327008365232    0.000730340585354261
        3       0.986493966243405    0.000184871077568705
        4       0.993202483991711    4.65199555944595e-05
        5       0.996589832769599    1.16688756544026e-05
        6       0.998292027147433    2.92215229613646e-06
        7       0.999145286557988    7.31159377154178e-07
        8          0.999572460931    1.82867799680686e-07
        9       0.999786184803889    4.57267127496053e-08
       10       0.999893080977268    1.14328996270174e-08
       11       0.999946537631102    2.85837764568697e-09
       12        0.99997326810061    7.14613479502191e-10
       13       0.999986633870396    1.78655867877353e-10
       14       0.999993316892757    4.46641612583676e-11
       15        0.99999665843075      1.116617909247e-11
       16       0.999998329221071    2.79143375081503e-12
       17       0.999999164587949    6.97886193279373e-13
       18       0.999999582277231    1.74527059471075e-13
       19       0.999999791180077    4.36317648677687e-14
       20       0.999999895652271    1.08801856413265e-14
       21       0.999999947786534    2.66453525910038e-15
       22       0.999999973302319     7.7715611723761e-16
       23       0.999999987857067    1.11022302462516e-16
       24       0.999999992428546                       0
牛顿重根情形迭代法
    IterNums     ApproximateSol         ErrorPrecision   
    ________    _________________    ____________________
       1         1.00384144150858    1.46999496749567e-05
       2         1.00000745433781     5.5566884427094e-11
       3        0.999999999988474                       0
复制成功

例3：求非线性方程在x_0%3D1附近的近似解：

f(x)%3D2e%5E%7B-x%7Dsinx%3D0

图4 各牛顿迭代法求解非线性方程误差精度下降曲线

其中下山法由于在第二次迭代时引入了下山因子0.25，故其下降曲线相对于牛顿法而言，保持了下降收敛的性质。

代码块
JavaScript
自动换行
复制代码
>> sol_downhill
sol_downhill = 
  包含以下字段的 struct:

           ApproximateSolution: 3.141592653589793
                ErrorPrecision: 1.058435733606881e-17
             NumerOfIterations: 7
                  NewtonMethod: "downhill"
    IntrationProcessInfomation: [7×3 double]
                   convergence: "收敛到满足精度的近似解"
                DownhillLambda: [1 0.250000000000000 1 1 1 1 1]
复制成功

科学计算 牛顿法 非线性方程求根 牛顿下山法 牛顿法重根情形 MATLAB程序设计
本文为我原创，未经授权禁止转载
cv16390567