数值分析之不动点迭代求多解（结合求多个近似解位置）

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

一、何为不动点迭代

不动点迭代和二分法的本质和区别在上一篇blog的开头就和大家说明了，作为简单迭代的一种，不动点迭代的关键在于寻找一个不定点（通过构建x=g(x)得到）。不动点迭代，只能对给定的一个近似解迭代找到解，所以需要结合求解近似值位置的算法，来找到区间上多个解的位置，然后依次对多个近似根迭代求值。

1.1 不动点迭代思想

这里只给出基本思想，不给出具体数学原理，有兴趣的可以自己搜其他的blog，讲的也很详细。
不动点迭代简而言之，是通过一个初始的近似根p，再构建一个x=g(x)，从p开始，利用公式Pn+1 = g(Pn)，逐次迭代逼近根，当这个根满足一个具体精度时，就认为求得的就是精确解。

matlab版算法：

1.2 求近似根位置

求解区间上所有的近似根，针对的是原函数f(x)=x**5-2*x-1。第一步就是将区间等分成若干段，这题中，为了满足两个样例，应该等分成8段，就是9个点。然后看是否满足以下2个条件任意之一，满足就认为这两点间存在解，近似根设为两点横坐标之和的一半。
（1）f(k-1)f(k)<0。#两个分别在x轴上下方，中间有根
（2）|y(k)|<eps_y and (y(k)-y(k-1))(f(k+1)-y(k))<0
#点靠x轴足够近，且曲线在点附近改变了斜率。

注意：由于截断误差和计算的误差，当y=f(x)非常陡峭的接近x轴时，则求解被认为是良性的，否则就是病态的。如果在p处有多个解，则只找到一个近似解。

matlab版算法：

二、题目及实现代码

2.1 题目

在[a,b]区间内寻找方程x**5-2*x-1=0的根的初始近似值位置，确定不动点迭代的初始点（可能有多个），然后使用不动点迭代法求方程的根（可能有多个根）。前后两次迭代的差的绝对值小于delta后停止迭代。求近似值位置的程序部分，epsilon取0.01。

2.2 输入输出格式

【输入形式】在屏幕上输入3个数，依次为区间左端点值a、右端点值b和所求根的精度值。各数间都以一个空格分隔。根据输入的所求根的精度值可求得delta.

【输出形式】每一行输出一个根（精确到小数点后3位）

2.3 样例

输入(1)

-1.2 1.5 3

输出(1)

-1.000
-0.519
1.291

输入(2)

-1.2 1.5 1

输出(2)

-1.012
-0.516
1.296

2.4 思路和要点

思路：将不动点迭代和求解近似根分成两个函数，先求解近似方根，将可能的根放入一个数列，然后带入不动点迭代的函数中，一次迭代求值。
求解近似根位置的函数：approot(arr_X) #arr_X是等分后区间的点集合
不动点迭代的函数： fixpt(Rlist, delta) #Rlist就是求出的近似根位置集合,delta就是输入的第三个数，精度要求。

要点：函数x**5-2*x-1=0在区间[-1.2,1.5]上有3个根，但是因为在三个点中，有的点的|g’§|>1,则Pn+1 = g(Pn)产生的序列对P发散，而只有满足收敛性才能迭代求正确解，所以要构建2个x=g(x),来满足3个点都要满足收敛性。具体图解如下。

由图可知，g1(x)在第二个点不满足收敛性，导数明显大于1，发散。所以在对第二个点迭代时，可以使用g2(x)。

2.5 代码

import numpy as np
def fun(x):
    return x**5-2*x-1
# 不动点迭代函数
def g_1(x):
    return np.sign(2*x+1)*pow(abs(2*x+1),1/5)
def g_2(x):
    return (x**5-1)/2

# 求解根的近似值位置
def approot(arr_X):
    np.seterr(invalid='ignore')
    eps = 0.01
    Rlist = []  # 存储根的列表
    arr_Y = fun(arr_X)  # 函数值
    yrange = max(arr_Y) - min(arr_Y)
    eps2 = yrange * eps
    n = len(arr_X)
    for k in range(1, n ):
        if arr_Y[k]*arr_Y[k-1] <= 0:
            Rlist.append((arr_X[k-1] + arr_X[k])/2)
        if k<n-3:
            #因为这里的(k+1)会溢出，导致for只执行到7，所以这里if条件排除溢出的这种情况
            s = (arr_Y[k] - arr_Y[k - 1]) * (arr_Y[k + 1] - arr_Y[k])
        else:
            s = 0
        if abs(arr_Y[k]) <= eps2 and s <= 0:
            Rlist.append(arr_X[k])
    return Rlist
 # 不动点迭代
def fixpt(Rlist, tol):
    np.seterr(invalid='ignore') #避免列表双标量无效时溢出的报错
    epsilon = np.finfo(np.float32).eps # machine epsilon
    outputlist = []
    k=0
    for num in Rlist:
        listpk = []
        listpk.append(num)
        while True:
            if k ==1:
                listpk.append(g_2(listpk[-1]))
            else:
                listpk.append(g_1(listpk[-1]))
            err = abs(listpk[-1] - listpk[-2])
            relerr = err / (abs(listpk[-1]) + epsilon)
            if err < tol or relerr < tol:
                break
        k += 1
        outputlist.append(listpk[-1])
    return outputlist

def main():
    n = np.array(input().split(),dtype = np.float32)
    a = n[0]
    b = n[1]
    delta = 10**(-int(n[2]))
    x = np.linspace(a, b, 9 ) #8个区间，有9个点
    Rlist = approot(x)
    root = fixpt(Rlist, delta)
    for i in root:
        print("%.3f" % i)

if __name__ == '__main__':
    main()

2.6 结果

样例1：

样例2：