数值分析之二分法、试值法 python

@ 数值分析之非线性方程求解

二分法、试值法的本质

二分法，试值法主要依靠在确定区间[a,b]上，f(a)f(b)<0来迭代找根，这个区间内只能有单解，然后缩小区间，逼近精确解，这统称为全局收敛法。
但是如果区间内有多个根，就要用不同区间来寻找根。而牛顿-拉夫森法和割线法能解决多解问题。
这类方法要求给定一个接近根的值保证收敛性，此称为局部收敛法，局部收敛速度大于全局收敛速度。
一些混合方法先采用全局收敛法，当迭代逼近根后切换局部收敛法。不动点迭代属于简单迭代，也只能求解单解问题，但是如果和求近似根的算法结合，就能处理多解问题，而这一综合算法的不动点迭代法将在下一篇展示。

(1) 二分法求利率

题目

如果在240个月内每月付款300美元，使用二分法在利率区间[a,b] 内，
求能够满足在这240个月之后,本金和利息总值达到50万美元的利率值，精确到小数点后第d位。

输入输出格式

【输入形式】在屏幕上输入3个数，依次为利率区间左端点值a、右端点值b和精确到小数点后d位。
各数间都以一个空格分隔。测试用例的输入满足：b>a>0, 1<=d<=8, d为正整数。

【输出形式】输出两行数据，第一行为迭代次数，第二行为求得的利率，保留d位小数。

举例

输入：

0.15 0.16 8

输出：

20
0.15753931

思路和要点

（1）精确到n位小数，需要一个相对误差的概念，一般来说，精确到小数点后n位，即delta=10e-n，要求最后的左右区间之差小于delta，也是迭代判定条件
（2）需要构建非线性方程，这里就是需要找到一个n（年利率）能带入方程与500000的差越小越好。这就需要知道本金利息是怎么计算的计算，公式如下。其中I是年利率，所以分摊到每月要除以12。

代码

// 二分法求利率
import numpy as np
def func(n):
    return 500000-300*12*(pow(1+n/12,240)-1)/n
    #这里是设定函数 f(n)是求最后得到的本兮和500000的差值
    #n/12的原因是n是年利率，分摊到每个月，要除以12 
def main():
    n = np.array(input().split(),dtype = np.float64)
    left = n[0]
    right = n[1]
    delta = 0.5*pow(10,-int(n[2]))
    count = 0   #迭代次数
    mid = (left+right)/2
    while (right-left)>delta: #非线性方程求近似解判定条件
        count += 1
        if func(mid)>0:
            left = mid
        elif func(mid)<0:
            right = mid
        elif func(mid)==0:  #完美解
            break
        mid = (left+right)/2
    print(count)
    print(round(mid,int(n[2])))
     
if __name__ =='__main__':
    main()

结果

(2)试值法法求利率

题目

如果在240个月内每月付款300美元，使用二分法在利率区间[a,b] 内，
求能够满足在这240个月之后,本金和利息总值达到50万美元的利率值，精确到小数点后第d位。

输入输出格式

【输入形式】在屏幕上输入3个数，依次为利率区间左端点值a、右端点值b和精确到小数点后d位。
各数间都以一个空格分隔。测试用例的输入满足：b>a>0, 1<=d<=8, d为正整数。

【输出形式】输出两行数据，第一行为迭代次数，第二行为求得的利率，保留d位小数。

举例

输入：

0.15 0.16 8

输出：

5
0.15753931

思路和要点

(1)试值法和二分法的共同点：试值法又叫试位法，和二分法的本质一样，都是依靠f(a)*f(b)<0来确定可导区间内至少存在一个值，通过缩小区间来求解，但是二分法收敛速度慢。所以试值法是对二分法的改进。
（2）改进之处：二分法是使用区间[a,b]的中点来迭代，试值法是寻找经过（a,f(a)）和（b,f(b)）的割线L与x轴的交点（c,0），则能得到一个更好的近似值，代码段中的fun_c(a,b)函数就是通过斜率来求c的值。

代码

// 试值法求利率
import numpy as np

def func(x):#带入x后的函数值
    return 500000-300 * 12 * (pow((1 + x / 12), 240) - 1) / x
def func_c(a, b):
#a为函数区间左端点，b是右端点
    return b - (func(b) * (b - a)) / (func(b) - func(a))
    #返回a,b两点连线与x轴的交点(c,0)的横坐标

def main():
    n = np.array(input().split(), dtype=np.float64)
    a = n[0]
    b = n[1]
    delta = 0.5 * 10 ** (-int(n[2]))
    #精确到小数点后n位：近似值的绝对误差小于小数点后n位上的半个单位
    epsilon = 0.0001
    count = 0
    if func(a) * func(b) > 0:#说明这段区间上没有解或有多个解
        return -1
    else:
        while True:
            count += 1  #count迭代次数
            c = func_c(a, b)
            if func(c)==0:  #c刚好就是所求解
                break
            elif func(a) * func(c) > 0:#解在c,b之间
                a = c
                c = func_c(a, b)
            elif func(a)*func(c)<0: #解在a,c之前
                b = c
                c = func_c(a, b)
            if abs(b - a) < delta: #区间范围<精度要求
                break
            if abs(func(c)) < epsilon: #函数值小于误差范围
                break       
    print(count)
    print(round(c, int(n[2])))

if __name__ == '__main__':
    main()

结果