数值分析之不动点迭代求多解（结合求多个近似解位置）
@数值分析之非线性方程求解
不动点迭代和二分法的本质和区别在上一篇blog的开头就和大家说明了，作为简单迭代的一种，不动点迭代的关键在于寻找一个不定点（通过构建x=g(x)得到）。不动点迭代，只能对给定的一个近似解迭代找到解，所以需要结合求解近似值位置的算法，来找到区间上多个解的位置，然后依次对多个近似根迭代求值。
1.1 不动点迭代思想这里只给出基本思想，不给出具体数学原理，有兴趣的可以自己搜其他的blog，讲的也很详细。
不动点迭代简而言之，是通过一个初始的近似根p，再构建一个x=g(x)，从p开始，利用公式Pn+1 = g(Pn)，逐次迭代逼近根，当这个根满足一个具体精度时，就认为求得的就是精确解。
matlab版算法：
求解区间上所有的近似根，针对的是原函数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版算法：
在[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
-1.2 1.5 1
输出(2)-1.012
-0.516
1.296
思路：将不动点迭代和求解近似根分成两个函数，先求解近似方根，将可能的根放入一个数列，然后带入不动点迭代的函数中，一次迭代求值。
求解近似根位置的函数：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)。
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：
作者：你要啥自行车
