问题描述

在[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

def f(x):
    return x ** 5 - 2 * x - 1

def g1(x):
    return (x**5-1)/2

def g2(x):
    result = (abs(2 * x + 1))**(1 / 5)
    if (2 * x - 1) < 0:
        return -result
    return result

def getEpsilon(x, epsilon):
    maxY = minY = x[0]
    for each in x:
        maxY = max(f(each), maxY)
        minY = min((f(each), minY))
    epsilon = (maxY - minY) * epsilon
    return epsilon

def getInitialVal(x, N, step, epsilon):
    initalVal = []
    for i in range(N + 1):
        y1, y2, y3 = f(x - step), f(x), f(x + step)
        if (y1 * y2 < 0) and (i != 0):
            initalVal.append((x + x - step) / 2)
        if ((y2 - y1) * (y3 - y2) < 0) and (abs(y2) < epsilon):
            initalVal.append(x)
        x += step

    return initalVal

def findFixedPoint(initalVal, delta,epsilon):
    points = []
    for each in initalVal:
        if (abs(g1(each)) < 1):
            points.append(iteration(each, g1, delta,epsilon))
        else:
            points.append(iteration(each, g2, delta,epsilon))
    return points

def iteration(p1, g, delta,epsilon):
    while True:
        p2 = g(p1)
        err =abs(p2-p1)
        relerr = err/(abs(p2)+epsilon)
        if err<delta or relerr<delta:
            return p2
        p1 = p2
                    
def main():
    a, b, c = input().split(' ')
    a = float(a)
    b = float(b)
    c = int(c)
    delta = 10 ** (-c)
    N = 8
    epsilon = 0.01
    step = (b - a) / N
    x = np.arange(a, b + epsilon, epsilon)
    
    epsilon2 = getEpsilon(x,epsilon)
    initalVal = getInitialVal(a, N, step, epsilon2)
    ans = findFixedPoint(initalVal, delta,epsilon)

    for each in ans:
        print('%.3f' % each)
        
if __name__ == '__main__':
    main()

迭代法的作用

许多复杂的求解问题，都可以转换成方程f(x)=0的求解问题。这一系列的解叫做方程的根。对于非线性方程的求解，在自变量范围内往往有多个解，我们将此变化区域分为多个小的子区间，对每个区间进行分别求解。我们在求解过程中，选取一个近似值或者近似区间，然后运用迭代方法逐步逼近真实解。
方程求根的常用迭代法有：二分法、不动点迭代、牛顿法、弦截法。

不动点迭代法

简单迭代法或基本迭代法又称不动点迭代法
1、不动点（FixedPoint）
首先来看一下什么是不动点：

换句话说，函数φ的不动点是y=φ(x)与y=x的交点，下图画出了函数y=cos(x)与y=x在区间[0,π/2]的交点，即cos(x)的不动点：

2、不动点迭代（Fixed Point Iteration）
不动点迭代的基本思想：

也就是说，为了求解方程f(x)=0，首先将方程转换为x=g(x)，然后初始化x0，循环迭代xi+1=g(xi)，直到满足收敛收件。
这里将方程f(x)=0转换为x=g(x)是很容易的，比如对于f(x)=x-cos(x)，求解f(x)=0即为求解x-cos(x)=0，即x=cos(x)，因此g(x)=cos(x)。
再例如对于方程：

可以等价为

还可以等价为

也就是说，将方程f(x)=0转换为x=g(x)有不同的方式，因此对方程f(x)=0来说，g(x)也不是唯一的。
3、不动点迭代的收敛性
这个迭代过程是很简单的，但这里有个关键性的问题：迭代收敛么？即经过N次迭代后是否会收敛于不动点？

通俗点讲，若要使不动点迭代收敛，则要求φ(x)在区间[a,b]上的函数值也在此区间内，另外还要求φ(x)的斜率绝对值不大于1。其证明过程比较复杂，有兴趣的可以查阅一些相关文献。

例题

求方程式：x³ - 0.165 × x² + 3.993 × 10^-4 = 0在（0，0.11）上的根

先看看不用迭代法计算的结果

from sympy import *
from sympy.abc import x

def func(x):
    return x**3 - 0.165*x**2 + 3.993*10**(-4)
result = solveset(func(x), x, Interval(0, 0.11))
print(result)

结果：

FiniteSet(0.0623775815137495)

约定一个误差，当误差小于某个数值的时候，迭代停止

代码：

xl = 0  #区间下限
xu = 0.11  #区间上限
x = (xl+xu)/2  #迭代初始值
x_list = [x]
i = 0

while True:
    x = x ** 3 - 0.165 * x ** 2 + 3.993 * 10 ** (-4) + x
    x_list.append(x)
    if len(x_list) > 1:
        i += 1
        error = abs((x_list[-1] - x_list[-2]) / x_list[-1])
        if error < 10**(-6):
            print(f'迭代第{i}次后，误差小于10^-6')
            break
    else:
        pass

结果：

迭代第777次后，误差小于10^-6
所求方程式的根为0.062370654088088

迭代至电脑默认误差为0

xl = 0  #区间下限
xu = 0.11  #区间上限
x = (xl+xu)/2  #迭代初始值
x_list = [x]
i = 0

while True:
    x = x ** 3 - 0.165 * x ** 2 + 3.993 * 10 ** (-4) + x
    x_list.append(x)
    if len(x_list) > 1:
        i += 1
        error = abs((x_list[-1] - x_list[-2]) / x_list[-1])
        if error == 0:
            print(f'迭代第{i}次后，误差为0')
            break
    else:
        pass

print(f'所求方程式的根为{x_list[-1]}')

结果：

迭代第3402次后，误差为0
所求方程式的根为0.062377581513749114

画迭代法

import matplotlib.pyplot as plt
xl = 0  #区间下限
xu = 0.11  #区间上限
x = (xl+xu)/2  #迭代初始值
x_list = [x]
i = 0

x_values = []
y_values = []
while True:
    x = x ** 3 - 0.165 * x ** 2 + 3.993 * 10 ** (-4) + x
    x_list.append(x)
    if len(x_list) > 1:
        i += 1
        error = abs((x_list[-1] - x_list[-2]) / x_list[-1])
        x_values.append(i)
        y_values.append(error)
        if error == 0:
            print(f'迭代第{i}次后，误差为0')
            break
    else:
        pass

print(f'所求方程式的根为{x_list[-1]}')

#设置绘图风格
plt.style.use('ggplot')
#处理中文乱码
plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']
#坐标轴负号的处理
plt.rcParams['axes.unicode_minus']=False
#横坐标是迭代次数
#纵坐标是误差值
plt.plot(x_values,
         y_values,
         color = 'steelblue', # 折线颜色
         marker = 'o', # 折线图中添加圆点
         markersize = 1, # 点的大小
         )
# 修改x轴和y轴标签
plt.xlabel('迭代次数')
plt.ylabel('误差值')
# 显示图形
plt.show()

结果：

迭代第3402次后，误差为0
所求方程式的根为0.062377581513749114

对比牛顿迭代法，牛顿迭代法要快很多，而且准确率也高
牛顿迭代法（Newton’s Method）迭代求根的Python程序

MATLAB基础

feval函数

用于求函数值

基本使用格式：y=feval(fhandle, x)    
%fhandle——函数表达式，x——变量值[y1, y2, ...] = feval(fhandle, x1,..., xn)

format long

设置输出格式，表示高精度输出

syms x

定义符号变量，用于占位符占位！

简单迭代法【不动点迭代】

大致讲解：

不动点可看成 $y=\phi (x)与y=x$ 的交点

MATLAB实现：

function [x,count]=Simple_iterative_method(x0,tol)
format long    %High precision output
error=1;    %The value of the initialization error
count=0;    %Initialization counter
while error>tol  %Shutdown criteria
    x=f(x0);   %Function to be solved
    error=abs(x-x0);  %Update the error of each step
    x0=x;  %Update iteration point
    count=count+1;   
end
end

%%Definition of function to be solved
function f=f(x)
f=(10/(4+x))^(1/2);
end

测试文件：

clear;  %Empty variable
clc;   %Clear screen
tol=1e-9;  %Set the required precision
x0=1.5;  %Exact solution of the point near
[x,count]=Simple_iterative_method(x0,tol);  %Use simple iterative algorithm
disp('The exact solution is as follows:')
disp(x);  
disp('Number of iterations required:')
disp(count);  

结果展示：

Newton 迭代法

matlab算法实现

function [x,count]=Newton_it(x0,tol)
format long
error=1;
count=0;
while error>tol
    x=x0-f(x0)/df(x0);  %Newton iteration scheme
    error=abs(x-x0);   
    x0=x;
    count=count+1;
end
end

%%Enter the function test below
function f=f(x0)
syms x;                %Defining symbolic variables
m(x)=x^2-115;
f=vpa(m(x0));      %Calculate the value of the symbolic variable
end

function df=df(x0)  %Calculating the first derivative of a function
syms x;                %Defining symbolic variables
m(x)=x^2-115;
df(x)=diff(m(x));
df=vpa(df(x0));
end

测试文件

clear;  %Empty variable
clc;   %Clear screen
x0=10;
tol=1e-6;
[x,count]=Newton_it(x0,tol);
disp('The exact solution is as follows:')
disp(x);
disp('Number of iterations required:')
disp(count);

结果展示：

作业

1. 用两种迭代方式求解一个多项式的根，多项式满足下面两个条件：
   （1）7次多项式
   （2）系数在(0,7] 之间
   要求：写出收敛阶、收敛速度、初值
2. 求出下面方程的根并作图$$x^2+{\left(\frac{5}{4}y-\sqrt{|x|}\right)}^2-4=0$$
   初值（这个方程的初值不好找！）

尝试解答：

1.如何生成一个随机整数向量（满足第二个条件叭）

round(rand(1,5)*6+1)

2.构造随机多项式

function f=ex6_1(x0,p0)     %Constructing polynomials
syms x;                %Defining symbolic variables
y(x)=poly2sym(p0);
disp(y(x));
f=vpa(y(x0));      %Calculate the value of the symbolic variable
end

p0=round(rand(1,8)*6+1);
%If you put it in a function, every time you call the function, a polynomial will be generated randomly
x0=2;
f=ex6_1(x0,p0);
disp(f);
x0=0;
f=ex6_1(x0,p0);
disp(f);

输出结果：

3.尝试使用简单迭代法

function [x,count]=Simple_iterative_method(x0,tol)
format long    %High precision output
error=1;    %The value of the initialization error
count=0;    %Initialization counter
while error>tol  %Shutdown criteria
    x=f(x0);   %Function to be solved
    error=abs(x-x0);  %Update the error of each step
    x0=x;  %Update iteration point
    count=count+1;   
end
end

%%Definition of function to be solved
function f=f(x0)
syms x;                %Defining symbolic variables
p0=round(rand(1,8)*6+1);  %Each time the function is called, a random polynomial is constructed
y(x)=poly2sym(p0);
f=vpa(y(x0));      %Calculate the value of the symbolic variable
end

clear;  %Empty variable
clc;   %Clear screen
tol=1e-9;  %Set the required precision
x0=0;  %Exact solution of the point near
[x,count]=Simple_iterative_method(x0,tol);  %Use simple iterative algorithm
disp('The exact solution is as follows:')
disp(x);  
disp('Number of iterations required:')
disp(count);  

输出结果（这是什么东西…）：

构造多项式并改写成 $x=\varphi (x)$ 的形式

function [f,df]=ex6_1(x0,p0)
syms x;                %Defining symbolic variablesS
y(x)=poly2sym(p0);
% disp('This is a randomly generated equation of degree 7')
% disp(y(x));
q(x)=((y(x)-p0(6)*x^2)/(-p0(6)))^(1/2);
% disp('This is phi(x)');
% disp(q(x));
f=vpa(q(x0));      %Calculate the value of the symbolic variable
df(x)=diff(q(x));
df=vpa(df(x0));
end

function [x,count]=Newton_it(x0,tol)
format long
error=1;
count=0;
p0=round(rand(1,8)*6+1);
while error>tol && count <50
    [f,df]=ex6_1(x0,p0);
    x=x0-f/df;  %Newton iteration scheme
    error=abs(x-x0);
    disp('error:')
    disp(error);
    x0=x;
    count=count+1;
end
end

clear;  %Empty variable
clc;   %Clear screen
x0=10;
tol=1e-3;
[x,count]=Newton_it(x0,tol);
disp('The exact solution is as follows:')
disp(x);
disp('Number of iterations required:')
disp(count);

function [x,count]=Simple_iterative_method(x0,tol)
format long    %High precision output
error=1;    %The value of the initialization error
count=0;    %Initialization counter
p0=round(rand(1,8)*6+1);
disp(p0);
while error>tol && count<20 %Shutdown criteria
    x=ex6_1(x0,p0);   %Function to be solved
    disp(x0);
    count=count+1;
    error=abs(x-x0);
    disp('it is error');
    disp(error);
    x0=x;
end
end

clear;  %Empty variable
clc;   %Clear screen
tol=1e-5;  %Set the required precision
x0=0;  %Exact solution of the point near
[x,count]=Simple_iterative_method(x0,tol);  %Use simple iterative algorithm
disp('The exact solution is as follows:')
disp(x);  
disp('Number of iterations required:')
disp(count);  

close all;
clear all;
clc;
ezplot('x^2+((5/4)*y-sqrt(abs(x)))^2-4=0',[-3,3])
grid on