非线性方程的数值解法-二分法,牛顿法,弦截法

1.二分法

1.1 定义

$考察有根区间[a,b],取中点x_0=(a+b)/2将它分成两半,$$假设中点x_0不是f(x)的零点,然后进行根的搜索,即检查$
$f(x_0)与f(a)是否同号,如果同号，说明所求的根x^{\ast }在x_0的右侧，这时令a_1=x_0,b_1=b;否则x^{\ast }必在x_0的左侧,$
$这时令a_1=a,b_1=x_0,对新的有根区间[a_1,b_1]施行同种过程,每次二分后，设取有根区间[a_k,b_k]的中点$
$$x_k=(a_k+b_k)/2$$
$作为根的近似值,则在二分过程中可以获得一个近似根序列$
$$x_0,x_1,x_2,\dots ,x_k,\dots ,$$
$该序列必以根x^{\ast }为极限,该方法称为二分法.$

1.2 Python实现二分法

from sympy.abc import x
from sympy import init_printing, evalf, exp, solve, Eq

init_printing(use_unicode=True)


def bisection(f_x, interval_start: (float, int), interval_end: (float, int), epsilon=5e-3):
    """
    二分法求解单变量非线性方程组f(x)=0在区间[a,b]上的准确到小数点后指定精度的近似解
    :param f_x: 单变量非线性函数f(x)
    :param interval_start: 求根区间始点
    :param interval_end: 求根区间终点
    :param epsilon: 小数点精度
    :return: 近似解x,二分迭代次数k
    """
    if f_x.subs(x, interval_start) * f_x.subs(x, interval_end) < 0:
        # 初始化二分次数，区间长度
        bisection_count = 0
        interval_length = abs(interval_start - interval_end)
        while 1:
            middle_point = (interval_start + interval_end) / 2
            bisection_count += 1
            value = f_x.subs(x, middle_point)
            if interval_length <= epsilon or value == 0:  # 达到预定精度epsilon或者等于精确解时结束循环
                return middle_point, bisection_count
            if value * f_x.subs(x, interval_start) < 0:
                interval_end = middle_point
            else:
                interval_start = middle_point
            interval_length = abs(interval_start - interval_end)
    else:
        raise Exception("function has no solution in the specified interval")

2.牛顿法

2.1 定义

$设已知方程f(x)=0有近似根x_k(假定f^{\prime }(x_k)\ne 0),将函数f(x)在点x_k处展开,有$
$$f(x)\approx f(x_k)+f^{\prime }(x_k)(x-x_k),$$
$于是方程f(x)=0可近似地表示为$
$$f(x_k)+f^{\prime }(x_k)(x-x_k)=0.$$
$这是个线性方程,记其根为x_{k+1}，则x_{k+1}的计算公式为$
$$x_{k+1}=x_k-\frac{f(x_k)}{f^{\prime }(x_k)},k=0,1,\dots ,$$
这就是牛顿法

2.2 Python实现牛顿法

def newtown_tangent(f_x, x0, true_root, epsilon=5e-4):
    """
    牛顿法求解单变量非线性方程组f(x)=0在区间[a,b]上的准确到小数点后指定精度的近似解
    :param f_x: 单变量非线性函数f(x)
    :param x0: 迭代初始近似值
    :param true_root: 精确解
    :param epsilon: 小数点精度
    :return: 近似解x,二分迭代次数k
    """
    f = f_x.subs(x, x0)
    diff_f_x = f_x.diff()
    diff_f = diff_f_x.subs(x, x0)
    iteration_count = 0
    while 1:
        x_k = x0 - f / diff_f
        iteration_count += 1
        if abs(x_k - true_root) < epsilon:
            return x_k, iteration_count
        x0 = x_k
        f = f_x.subs(x, x0)
        diff_f = diff_f_x.subs(x, x0)

3.弦截法

3.1 定义

$设x_k,x_{k+1}是方程f(x)=0的近似根，$$利用f(x_k)，f(x_{k-1})构造一次插值多项式p_1(x),并用p_1(x)=0$
$的根作为方程f(x)=0的新的近似根x_{k+1}.由于$
$$p_1(x)=f(x_k)+\frac{f(x_k)-f(x_{k-1})}{x_k-x_{k-1}}(x-x_k),(1)$$
$因此有$
$$x_{k+1}=x_k-\frac{f(x_k)}{f(x_k)-f(x_{k-1})}(x_k-x_{k-1}),k=0,1,\dots ,(2)$$
$这样导出的迭代公式(2)可以看做牛顿公式$
$$x_{k+1}=x_k-\frac{f(x_k)}{f^{\prime }(x_k)}$$
$中导数f^{\prime }(x_k)用差商\frac{f(x_k)-f(x_{k-1})}{x_k-x_{k-1}}取代的结果.$
$该方法称为弦截法.$

3.2 Python实现弦截法

def secant(f_x, x0, x1, true_root, epsilon=5e-4):
    """
    弦截法求解单变量非线性方程组f(x)=0在区间[a,b]上的准确到小数点后指定精度的近似解
    :param f_x: 单变量非线性函数f(x)
    :param x0: 第一个迭代初始近似值
    :param x1: 第二个迭代初始近似值
    :param true_root: 精确解
    :param epsilon: 小数点精度
    :return: 近似解x,二分迭代次数k
    """
    f_x0 = f_x.subs(x, x0)
    f_x1 = f_x.subs(x, x1)
    difference_quotient = (f_x1 - f_x0) / (x1 - x0)
    iteration_count = 0
    while 1:
        x_k = x1 - f_x1 / difference_quotient
        iteration_count += 1
        if abs(x_k - true_root) < epsilon:
            return x_k, iteration_count
        x0 = x1
        x1 = x_k
        f_x0 = f_x.subs(x, x0)
        f_x1 = f_x.subs(x, x1)
        difference_quotient = (f_x1 - f_x0) / (x1 - x0)

4.测试

if __name__ == '__main__':
    # 第一组测试函数，来源详见李庆扬数值分析第5版P214，e.g.2
    f = x ** 3 - x - 1
    appropriate_root1, iteration_count1 = bisection(f, 1.0, 1.5)
    print("二分法求方程近似解x={},迭代次数={}".format(appropriate_root1,iteration_count1))
    # 二分法测试区间中点恰好是精确解时的情况
    f = x - 1
    appropriate_root2, iteration_count2 = bisection(f, 0.0, 2.0)
    print("测试区间中点恰好是精确解时,近似解x={},迭代次数={}".format(appropriate_root2, iteration_count2))
    # 第二组测试函数，来源详见李庆扬数值分析第5版P238，e.g.7
    f = x ** 3 - 3 * x - 1
    initial_value = 2
    true_solution = 1.87938524
    root1, count1 = newtown_tangent(f, initial_value, true_solution)
    print("牛顿切线法近似解x={},迭代次数={}".format(root1.evalf(), count1))
    initial_value0, initial_value1 = 2, 1.9
    root2, count2 = secant(f, initial_value0, initial_value1, true_solution)
    print("弦截法近似解x={},迭代次数={}".format(root2.evalf(), count2))
    # 第三组测试函数，来源详见李庆扬数值分析第5版P224-225
    c = 115
    f = x ** 2 - c
    true_solution = pow(c, 0.5)
    initial_value = 11
    solution, count = newtown_tangent(f, initial_value, true_solution, epsilon=1e-6)
    print("牛顿切线法近似解x={},迭代次数={}".format(solution.evalf(), count))
    # 第四组测试函数，来源详见李庆扬数值分析第5版P229 e.g.10
    f = x * exp(x) - 1
    root = solve(Eq(f, 0), x)[0]
    print("调用sympy函数库解x={}".format(root))
    initial_value0, initial_value1 = 0.5, 0.6
    root, count = secant(f, initial_value0, initial_value1, root)
    print("牛顿切线法近似解x={},迭代次数={}".format(root.evalf(), count))

5.运行结果