非线性方程的解法1-3， 1-4

一. 不动点迭代法

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

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

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

【样例1输入】

-1.2 1.5 3

【样例1输出】

-1.000

-0.519

1.291

【样例1说明】输入：左端点a值为-1.2，右端点b值为1.5，前后两次迭代的差的绝对值小于delta=10**(-3)后停止迭代。输出：从小到大顺序输出三个根的值。

【评分标准】根据输入得到的输出准确

import numpy as np
# 求x的函数值，函数为f（x） = x**5 -2*x -1
def f(x):
    y = x ** 5 - 2 * x - 1
    return y
# 确定根的初始近似值
def approot(x, epsilon):
    y = f(x)
    y_range =np.max(y) - np.min(y)
    epsilon2 = y_range * epsilon
    x = np.append(x, x[-1])
    y = np.append(y, y[-1])
    #初始化r为空的array
    r = np.array([], dtype=np.double)
    for k in range(1, x.size-1):
        if y[k-1]*y[k]<=0:
           r = np.append(r, (x[k-1]+x[k])/2)
        s =(y[k]-y[k-1])*(y[k+1]-y[k])
        if np.abs(y[k])<epsilon2 and (s<=0):
           r = np.append(r, x[k])
    return r

def g1(x):
    value = pow(abs(x * 2 + 1), 1 / 5)
    if x < 0:
       value = -value
    return value
def g2(x):
    value = (x ** 5 - 1) / 2
    return value

def fix_point(x0, delta, max_iter, k):
     epsilon = np.finfo(np.float32).eps
     #根序列
     r_seq = np.zeros(max_iter, dtype=np.double)
     r_seq[0] = x0
     #初始化
     r = 0
     n = 0
     err = 0
     for n in range(1, max_iter):
         if k == 0 or k == 2:
             r_seq[n] = g1(r_seq[n-1])
         else:
             r_seq[n] = g2(r_seq[n-1])
         err = np.abs(r_seq[n] - r_seq[n - 1])
         rel_err = err/np.abs(r_seq[n] + epsilon)
         r = r_seq[n]
         if(err < delta)or(rel_err < delta):
             break
     return r, n, err

def main():
        # 创建一个初始值的列表
        arr = input().split()
        arr = np.array(arr, dtype = float)
        #求根的近似值位置
        a = arr[0]
        b = arr[1]
        #在区间[a, b]内取多个点
        #np.linspace创建等差数列
        x = np.linspace(a, b, 9)
        #epsilon是公差
        epsilon = 0.01
        # 有几个根，就有几个近似位置
        r_app = approot(x, epsilon)
        #if r_app is not None: print(r_app)# 输出近似位置
        # 对初始值进行迭代求根
        delta = 0.5 * 10 ** (-int(arr[2]))#精度
        max_iter = 50 #最大迭代次数
        r = np.zeros(r_app.size, dtype=np.double)
        for k in range(r_app.size):
            r[k],_,_ = fix_point(r_app[k], delta, max_iter, k)
            print("%.3f" %r[k])
if __name__ == '__main__':
    main() 

二 . 牛顿法解投射问题

【题目简述】牛顿法解投射体问题（非线性方程求解）

【问题描述】在考虑空气阻力情况下，求解投射体撞击地面时经过的时间和水平飞行行程，其中：y=f(t)=9600*(1-e**(-t/15.0)) - 480t；x=r(t)=2400(1-e**(-t/15.0))。

【输入形式】在屏幕上输入3个数，分别表示起始值、前后两次迭代的差的绝对值精度和f(t)函数值的精度。各数间都以一个空格分隔。

【输出形式】第一行输出投射体撞击地面时经过的时间，第二行输出水平飞行行程，精确到小数点后5位。

【样例1输入】

8 1 1

【样例1输出】

9.08955

1090.69211

【样例1说明】输入：起始值为8、前后两次迭代的差的绝对值精度为0.1和f(t)函数值的精度为0.1。输出：第一行输出投射体撞击地面时经过的时间为9.08955秒，第二行输出水平飞行行程为1090.69211ft。

import numpy as np
import math
def f(t):
    #return 垂直行程
    #e = 2.71828
    y = 9600*(1-math.e**(-t/15.0)) - 480*t
    return y
def r(t):
    # return 水平行程
    #e = 2.71828
    x=2400*(1-math.e**(-t/15.0))
    return x
def ff(t):
    #return f(t)一阶导数的值
    yy = 640*math.e**(-t/15.0)-480
    return yy
def newton(p0, delta, eps, max_iter):
    for k in range(max_iter):
        p1 = p0 - f(p0)/ff(p0)
        err = np.abs(p1-p0)
        rel_err = 2*err/(np.abs(p1)+delta)
        p0 = p1
        if(err<delta or np.abs(f(p0)) <eps or rel_err < delta):
            break
    x= r(p0)
    return p0,x
    """
    y=1
    rel=1
    err =1
    while(err>=delta and abs(y) >=eps and rel >= delta ):
        p1 = p0-f(p0)/ff(p0)
        err = abs(p1 - p0);
        rel = 2*err/(abs(p1)+delta)
        y=f(p0)-f(p1)
        p0 = p1
        #y = f(p0)
        x= r(p0)
    return p0,x
    """
def main():
    p0, n, m = input().split()
    #起始值, 前后两次迭代的差的绝对值精度, f(t)函数值的精度
    p0= float(p0)
    n = int(n)
    m = int(m)
    delta=10 ** (-n) #前后两次迭代的差的绝对值精度
    eps=10 ** (-m) #f(t)函数值的精度
    max_iter = 50
    ans,x = newton(p0, delta, eps, max_iter)
    print("%.5f" % ans) #时间
    print("%.5f" % x) #水平距离
if __name__ == '__main__':
    main()