python 牛顿迭代法求解非线性方程有限区间内所有实数数值解

1. 引言

上一篇文章python 二分法求解非线性方程所有实数解，初步实现了对非线性方程解的查找。由于采用的是二分法，收敛速度较慢，这里改进为牛顿迭代，并对收敛速度进行对比。

2.程序

2.1 调用库

import numpy as np
import matplotlib.pyplot as plt

2.2 求解程序

def solveFun(fun, min_x=0, max_x=1, n_step=1000, max_lter=64, max_error= 1e-6, plot=False, print_log=False):
    """
    :param fun: 需要求解的函数
    :param min_x: 自变量求解区间下限
    :param max_x: 自变量求解区间上限
    :param n_step: 查找解的步数，区间长度/步长
    :param max_lter: 最大迭代次数
    :param max_error: 容许相对误差
    :param plot: True绘制曲线图， False不绘制曲线
    :param print_log: True输出求解信息，False不输出求解信息
    :return: numpy.array()对象，包含解集
    """
    x_list = np.linspace(min_x, max_x, n_step)
    y_list = np.zeros(n_step)
    solve_list = []
    N_solve=0
    for i in range(n_step):
        y_list[i] = fun(x_list[i])
        if  0 < i and y_list[i] * y_list[i - 1] <= 0:
            N_solve+= 1
            temp1_x = x_list[i]
            temp2_x = x_list[i - 1]
            temp1_y = y_list[i]
            temp2_y = y_list[i - 1]
            dy= (temp1_y- temp2_y)/(temp1_x- temp2_x)
            error_ = abs((temp1_x - temp2_x) / (abs(temp1_x) + 1e-10)) # 给分母加上一个较小数值，以免出现除零的情况
            lter_n = 0
            if print_log:
                print("找到第 %d个解位于 %11.4E 附近，进入迭代...\n"
                      "     迭代次数         相对误差                x               fun" % (N_solve, temp1_x))

            while lter_n < max_lter and error_ > max_error:
                lter_n += 1
                temp2_x= temp1_x
                temp1_x= temp2_x- temp1_y/dy
                temp2_y= temp1_y
                temp1_y = fun(temp1_x)
                dy= (temp1_y- temp2_y)/(temp1_x- temp2_x)
                error_ = abs((temp1_x - temp2_x) / (abs(temp1_x) + 1e-10))
                if print_log:
                    print(" %8d   %16.4E    %16.4E   %16.4E" % (lter_n, error_, temp1_x, temp1_y))
                temp2_x = temp1_x
            if lter_n < max_lter:
                solve_list.append(temp1_x)
                if plot:
                    plt.plot(temp1_x, temp1_y, marker=".", markersize=16, color="r")
                    # plt.text(temp1_x, temp1_y, "%6.4E" % temp1_x, color="k")
                # print("迭代次数: %2d; res %d: kr= %6.4E, form4=%6.4E, error= %6.4E" % (
                #     max_in, len(res_list), temp1_kr, temp1_f4, error_i))
            else:
                print("达到最大迭代次数，相对误差 %16.4E" % error_)
        # print(res_list[i])
    if plot:
        plt.plot(x_list, y_list, color="b")
        plt.show()
    return np.array(solve_list)

3. 测试

3.1 计算测试

测试程序：

def myFun(x):
    return np.sin(3 * x**0.7) + 0.5 * np.sin(2 * x + 1)


print(solveFun(fun=myFun, min_x=0, max_x=20.0, print_log= True, plot= True)) #

命令行输出：

找到第 1个解位于  1.0811E+00 附近，进入迭代...
     迭代次数         相对误差                x               fun
        1         1.1309E-02          1.0690E+00        -2.8511E-05
        2         8.7314E-06          1.0690E+00         3.5370E-08
        3         1.0819E-08          1.0690E+00        -3.2070E-14
找到第 2个解位于  2.8028E+00 附近，进入迭代...
     迭代次数         相对误差                x               fun
        1         6.8830E-03          2.7836E+00         3.5998E-06
        2         5.2038E-07          2.7836E+00        -6.0916E-09

... ...


找到第 7个解位于  1.6997E+01 附近，进入迭代...
     迭代次数         相对误差                x               fun
        1         5.8115E-04          1.6987E+01        -3.3137E-05
        2         1.0896E-06          1.6987E+01         6.1904E-08
        3         2.0317E-09          1.6987E+01        -2.1566E-13
[ 1.0689822   2.78364163  4.77207543  8.20017307 10.57158515 14.35831954
 16.98710642]

Process finished with exit code 0

绘图输出：

 

3.2 与二分法对比

二分法迭代输出：

找到第 7个解位于( 1.6977E+01,  1.6997E+01)，进入迭代...
     迭代次数         相对误差                x               fun
        1         8.8365E-04          1.6992E+01         2.1423E-04
        2         1.4730E-04          1.6989E+01        -8.7551E-03
        3         7.3654E-05          1.6988E+01        -4.2725E-03
        4         3.6828E-05          1.6988E+01        -2.0297E-03
        5         1.8415E-05          1.6987E+01        -9.0785E-04
        6         9.2073E-06          1.6987E+01        -3.4684E-04
        7         4.6037E-06          1.6987E+01        -6.6311E-05
        8         2.3018E-06          1.6987E+01         7.3960E-05
        9         1.1509E-06          1.6987E+01         3.8239E-06
       10         5.7546E-07          1.6987E+01        -3.1244E-05

牛顿迭代法输出：

找到第 7个解位于  1.6997E+01 附近，进入迭代...
     迭代次数         相对误差                x               fun
        1         5.8115E-04          1.6987E+01        -3.3137E-05
        2         1.0896E-06          1.6987E+01         6.1904E-08
        3         2.0317E-09          1.6987E+01        -2.1566E-13

可以观察到牛顿迭代法收敛速度远远优于二分法。二分法进行了10次迭代才将误差降低到5.75E-7， 而牛顿迭代仅用了三次迭代就将误差降低到2.03E-9。