二分法


```python
from sympy import *
import math
```


```python
def bisection(f, a, b, eps, eta=1e-16, verbose=False):
    """二分法求根
    
    对方程 f(x) = 0 在区间 [a, b]，使用二分法求根。
    做 ceil((log(((b - a) / eps), 2) - 1) 次迭代，使结果满足精度 eps。
    
    实现参考: 数值分析[谷根代，杨晓忠 等著]2011年版.P18.算法1

    Args:
        f: function, 一元函数，表示要求根的方程: f(x) = 0.
        a, b: float, 有根区间 [a, b] 的端点.
        eps: float, 给定精度.
        
        eta: float, 当 abs(f(x)) <= eta 时停止计算, default 1e-16.
        verbose: bool, 打印出二分法计算的表格, default False.

    Returns:
        (x_final, N)

        x_final: float, 二分法求得的近似根
        N: int, 迭代次数

    Raises:
        ValueError: 给定的 f(a) * f(b) < 0 时无法求解，抛出异常
    """
    if f(a) * f(b) > 0:
        raise ValueError("rootless interval: f(a) * f(b) < 0")

    N = math.ceil((log(((b - a) / eps), 2) - 1).evalf(5))

    if verbose:
        print(f'N = {N}\n')
    
    x = (a + b) / 2
    
    if verbose:
        print(f'n \t (a, b) \t f(x_n)')
        print('-'*35)
    
    for n in range(N+1):
        if abs(f(x)) <= eta:  # f(x) == 0
            break

        if verbose:
            print(f'{n} \t ({a}, {b}) \t f({x})={f(x)}')
        
        if f(x) * f(a) < 0:
            b = x
        else:
            a = x
        x = (a + b) / 2
        n += 1
        
    if verbose:
        print(f'{n} \t ({a}, {b}) \t -')
    
    x_final = (a + b) / 2
    if verbose:
        print(f'\nresult: x = ({a}+{b})/2 = {x_final}')

    return x_final, N+1
    
```


```python
def f(x):
    return 2 * exp(-x) - sin(x)

bisection(f, 0, 1, 0.0005, verbose=True)
```

    N = 10
    
    n 	 (a, b) 	 f(x_n)
    -----------------------------------
    0 	 (0, 1) 	 f(0.5)=0.733635780821064
    1 	 (0.5, 1) 	 f(0.75)=0.263094345458695
    2 	 (0.75, 1) 	 f(0.875)=0.0661805371209897
    3 	 (0.875, 1) 	 f(0.9375)=-0.0228698549070950
    4 	 (0.875, 0.9375) 	 f(0.90625)=0.0208763999947352
    5 	 (0.90625, 0.9375) 	 f(0.921875)=-0.00119109771321235
    6 	 (0.90625, 0.921875) 	 f(0.9140625)=0.00979401298210625
    7 	 (0.9140625, 0.921875) 	 f(0.91796875)=0.00428930379438952
    8 	 (0.91796875, 0.921875) 	 f(0.919921875)=0.00154606529293533
    9 	 (0.919921875, 0.921875) 	 f(0.9208984375)=0.000176724441991016
    10 	 (0.9208984375, 0.921875) 	 f(0.92138671875)=-0.000507376461447939
    11 	 (0.9208984375, 0.92138671875) 	 -
    
    result: x = (0.9208984375+0.92138671875)/2 = 0.921142578125





    (0.921142578125, 11)




```python
def dichotomy(f, a, b, eps, eta=1e-16, verbose=False):
    """二分法求根

    对方程 f(x) = 0 在区间 [a, b]，使用二分法求根。
    一直迭代到 abs(b - a) <= eps 为止。

    实现参考：《实验二  非线性方程求根》

    Args:
        f: function, 一元函数，表示要求根的方程: f(x) = 0.
        a, b: float, 有根区间 [a, b] 的端点.
        eps: float, 根的容许误差.
        
        eta: float, abs(f(x)) 的容许误差, default 1e-16.
        verbose: bool, 打印出二分法计算的表格, default False.

    Returns:
        (x_final, N)

        x_final: float, 二分法求得的近似根
        N: int, 迭代次数

    Raises:
        ValueError: 给定的 f(a) * f(b) < 0 时无法求解，抛出异常
    """
    if f(a) * f(b) > 0:
        raise ValueError("rootless interval: f(a) * f(b) < 0")
    
    if verbose:
        print(f'n \t (a, b) \t f(x_n)')
        print('-'*35)
    
    n = 0
    while abs(b - a) > eps:
        x = (a + b) / 2
        fx = f(x)
        if verbose:
            print(f"{n}\t {a, b}\t f({x})={fx}")

        if abs(fx) <= eta:
            break
        
        if f(a) * fx < 0:
            b = x
        else:
            a = x
        
        n += 1

    x_final = (a + b) / 2
    if verbose:
        print(f"{n}\t {a, b}\t -")
        print(f'\nresult: x = ({a}+{b})/2 = {x_final}')
    
    return x_final, n

```


```python
dichotomy(f, 0, 1, 0.0005, 0)
```




    (0.921142578125, 11)




```python
def fixed_point_iter(phi, x_0, max_steps=25, verbose=False):
    """不动点迭代

    Args: 
        phi: function, 迭代函数
        x_0: float, 初值
        
        max_steps: int, 最大迭代次数
        verbose: bool, 打印出每一步的值，default False.

    Returns: 
        x_final: float, 最终的近似根 x 
    """
    x = x_0
    for i in range(max_steps):
        if verbose:
            print(f'x_{i} \t {x}')
        x = phi(x)
    return x
```


```python
fixed_point_iter(lambda x: 1 + 1 / x**2, 1.5, max_steps=10, verbose=True)
```

    x_0 	 1.5
    x_1 	 1.4444444444444444
    x_2 	 1.4792899408284024
    x_3 	 1.456976
    x_4 	 1.4710805833200253
    x_5 	 1.4620905354712408
    x_6 	 1.4677905760195855
    x_7 	 1.464164380462178
    x_8 	 1.4664663557170745
    x_9 	 1.465003040566855





    1.4659324390818347




```python
def newton_iter(f, x_0, eps=0, eta=0, df=None, max_steps=20, frac=False, verbose=False):
    """牛顿迭代法

    Args: 
        f: function, 迭代函数
        x_0: float, 初值
        eps: float, 根的容许误差, default 0.
        eta: float, abs(f(x)) 的容许误差, default 0.
        
        df: function, f 的导函数, 默认 None 表示自动调用 sympy.diff 求导(会导致后续迭代中使用分数运算)。
        max_steps: int, 最大迭代次数, default 20.
        frac: bool, True 则输出分数(仅对 df=None 时生效)，否则使用 float, default False.
        verbose: bool, 打印出每一步的值, default False.

    Returns: 
        x_final: float, 最终的近似根 x
    """

    if df == None:
        __x = Symbol('x')
        __df = diff(f(__x), __x)
        df = lambda x: __df.subs(__x, x)

    x = x_0
    for i in range(max_steps):
        if verbose:
            print(i, x)
        
        x_next = x - f(x) / df(x)
        
        if abs(x_next - x) <= eps or abs(f(x)) <= eta:
           break
        
        x = x_next
    
    if not frac:
        x = float(x)
    return x
```


```python
def f(x):
    return x ** 3 - 2 * x - 5

def df(x):
    return 3 * x ** 2 - 2

# newton_iter(f, 2, eps=0.5e-3, frac=True, verbose=True)
newton_iter(f, 2, eps=0.5e-3, df=df, frac=True, verbose=True)
```

    0 2
    1 2.1
    2 2.094568121104185





    2.094568121104185




```python
def single_point_truncation(f, x_0, x_1, eps=0, eta=0, max_steps=20, verbose=False):
    """单点弦截法

    Args: 
        f: function, 迭代函数
        x_0, x_1: float, 初值
        eps: float, 根的容许误差, default 0.
        eta: float, abs(f(x)) 的容许误差, default 0.
        
        max_steps: int, 最大迭代次数, default 20.
        verbose: bool, 打印出每一步的值, default False.

    Returns: 
        x_final: float, 最终的近似根 x
    """

    f_x0 = f(x_0)
    
    x = x_1
    for i in range(max_steps):
        if verbose:
            print(i, x)
        
        x_next = x - f(x) / (f(x) - f_x0) * (x - x_0)
        
        if abs(x_next - x) <= eps or abs(f(x)) <= eta:
           break
        
        x = x_next
    
    return x
```


```python
def f(x):
    return x ** 3 - 2 * x - 5

single_point_truncation(f, 2, 1, eps=0.5e-3, verbose=True)
```

    0 1
    1 2.2
    2 2.088967971530249
    3 2.094861151990966





    2.094861151990966




```python
def secant_method(f, x0, x1, eps=0, eta=0, max_steps=20, verbose=False):
    """两点弦截法（割线法）

    Args: 
        f: function, 迭代函数
        x0, x1: float, 初值
        eps: float, 根的容许误差, default 0.
        eta: float, abs(f(x)) 的容许误差, default 0.
        
        max_steps: int, 最大迭代次数, default 20.
        verbose: bool, 打印出每一步的值, default False.

    Returns: 
        x_final: float, 最终的近似根 x
    """
    for i in range(max_steps):
        if verbose:
            print(i, x1)
        
        x2 = x1 - f(x1) * (x1 - x0) / (f(x1) - f(x0))
        x0, x1 = x1, x2

        if abs(x1 - x0) <= eps or abs(f(x1)) <= eta:
           break
            
    return x1
```


```python
secant_method(lambda x: x ** 2 - 612, 10, 30, max_steps=5, verbose=True) # Root: 24.738633748750722
```

    0 30
    1 22.8
    2 24.545454545454547
    3 24.746543778801843
    4 24.73860275369709





    24.738633748750722



题目 
求方程 $f(x)=x^3+x^2-3x-3=0$ 在 1.5 附近的根.（误差限为 $\epsilon=1e-6$, $\eta=1e-9$）



```python
def f(x):
    return x ** 3 + x ** 2 - 3 * x - 3

x0 = 1.5
x1 = 2
eps = 1e-6
eta = 1e-9

result = {}
print("牛顿迭代法:")
result["牛顿迭代法"] = newton_iter(f, x0, eps=eps, eta=eta, verbose=True)
print("单点弦截法:")
result["单点弦截法"] = single_point_truncation(f, x0, x1, eps=eps, eta=eta, verbose=True)
print("两点弦截法:")
result["两点弦截法"] = secant_method(f, x0, x1, eps=eps, eta=eta, verbose=True)

print("\nresult:")
for k in result:
    print(f'{k}: {result[k]}')
```

    牛顿迭代法:
    0 1.5
    1 1.77777777777778
    2 1.73336066694000
    3 1.73205192940947
    4 1.73205080756970
    单点弦截法:
    0 2
    1 1.6923076923076923
    2 1.7390156515180086
    3 1.7308625826467308
    4 1.7322544663053168
    5 1.7320159286895187
    6 1.7320567817872679
    7 1.7320497843005194
    8 1.732050982835706
    两点弦截法:
    0 2
    1 1.6923076923076923
    2 1.7257977285018928
    3 1.7322172842612025
    4 1.732050123979108
    
    result:
    牛顿迭代法: 1.7320508075697012
    单点弦截法: 1.732050982835706
    两点弦截法: 1.7320508074943775



```python

```

0ae96761a77ef9afdd12c0a5124b138c9ccb6fb3
ipynb
Jupyter Notebook
ex2/src/ex2.ipynb
cdfmlr/NumericalAnalysis
b752af60c1f3202d53cc3c76e66419a885250ed9
Apache-2.0
ex2/src/ex2.ipynb
cdfmlr/NumericalAnalysis
b752af60c1f3202d53cc3c76e66419a885250ed9
Apache-2.0
ex2/src/ex2.ipynb
cdfmlr/NumericalAnalysis
b752af60c1f3202d53cc3c76e66419a885250ed9
Apache-2.0
2021-08-15T01:34:35.000Z
2021-08-15T01:34:35.000Z
Qwen/Qwen-72B
1. YES
2. YES