非线性方程的数值解法

最新推荐文章于 2023-10-23 10:27:13 发布

21e23f234

最新推荐文章于 2023-10-23 10:27:13 发布

阅读量617

点赞数 1

文章标签：笔记算法 python

本文链接：https://blog.csdn.net/qq_28247201/article/details/132943885

版权

教材：数值计算方法第3版马东升
非线性方程：高次代数方程和超越方程

方程的根

定义 非线性方程 f(x)=0，存在x使f(x)=0 ，则称 x*为方程的根，又称为函数f(x) 的零点。
方程的根：实根、虚根。全局的根、局部的根。单根、重根。

单根

函数f(x)对于数x有 f(x)=0，但 f´(x*)≠0 则称x*为方程的单根。

重根

若f(x)满足f(x)=(x-x*)mβ(x)，其中β(x)≠0，则称x是方程f(x)=0的m重根和函数f(x)的m重零点。当m=1时，x是方程的单根和函数f(x)的单零点。

介值定理

若函数f(x)在[a,b] 连续，且f(a)*f(b)<0，则方程f(x)=0 在(a,b)内至少有一个实根。将[a,b]称为f(x)的有根区间。

一、初始近似值的搜索

1.逐步搜索法

假设f(x)在区间[a,b]内有一个实根x*,若b-a较小，则可在(a,b)上任取一点x0作为初始近似根。一般情形，可用逐步搜索法

给定有根区间 (a, b)，其中 f(a) < 0 且 f(b) > 0，以及步长 h。
初始化变量 x 为左端点 a，即 x = a。
开始迭代：

计算函数值 f(x)。
如果 f(x) 接近零（可以设置一个小的阈值来判断），则认为 x 是一个根，停止搜索。
如果 f(x) < 0，则说明根位于当前点的右侧，需要向右跨一步。更新 x = x + h。
如果 f(x) > 0，则说明根位于当前点的左侧，需要缩小有根区间。确定一个新的有根区间 (x - h, x)，其宽度等于步长 h。可以选择不同的策略来确定新的有根区间，例如二分法等。

重复步骤 3，直到找到满足条件的根或达到搜索的终止条件。

逐步搜索法确定步长h是个关键。只要步长h取得足够小，利用这种方法可以得到具有任意精度的近似根。不过当减小步长时，搜索的步数相应增多，从而使计算量加大。因此，逐步搜索法一般用于初步确定根的位置。

测试

f(x)=𝑥^3−4𝑥ln(𝑥+2)−1,在区间[0,4]中找零点

import math

# 定义要搜索的函数 f(x)


def f(x):
    return x**3-4*x*math.log(x+2)-1


# 初始化搜索参数
a = 0.0  # 区间左端点
b = 4.0  # 区间右端点

threshold = 0.01  # 阈值，用于判断 f(x) 是否接近零
max_iterations = 10000  # 最大迭代次数
h = (b-a)/max_iterations  # 步长

# 初始化迭代变量
x = a
iteration = 1

# 打印表头
print("| 迭代次数 | 当前x   | f(x)     | 区间 (a, b)   |")
print("|---------|--------|----------|---------------|")

with open("root-finding-results.txt", "a") as file:
    file.write("| 迭代次数 | 当前x   | f(x)     | 区间 (a, b)   |\n")
    file.write("|---------|--------|----------|---------------|\n")
    # 开始迭代搜索
    while iteration <= max_iterations:
        # 计算当前函数值
        fx = f(x)

        if fx < 0:
            # 根位于当前点的右侧，向右跨一步
            x = x + h
            a = x  # 更新左端点
        elif fx > 0:
            # 根位于当前点的左侧，缩小有根区间

            b = x  # 更新右端点
            a = x-h # 更新左端点
            if(abs(f(x)) <= threshold):
                print(f"| {iteration}\t| {x:.4f}\t| {fx:.4f}\t| ({a:.4f}, {b:.4f})\t|")
                file.write(f"| {iteration}\t| {x:.4f}\t| {fx:.4f}\t| ({a:.4f}, {b:.4f})\t|\n")
                break

        print(f"| {iteration}\t| {x:.4f}\t| {fx:.4f}\t| ({a:.4f}, {b:.4f})\t|")
        file.write(f"| {iteration}\t| {x:.4f}\t| {fx:.4f}\t| ({a:.4f}, {b:.4f})\t|\n")

        iteration += 1

    # 打印最终结果
    print("\n搜索结束:")
    if abs(f(x)) <= threshold:
        print(f"找到根：x = {x:.4f}, f(x) = {f(x):.4f}")
    else:
        print("未找到满足条件的根")

部分输出

迭代次数	当前x	f(x)	区间 (a, b)
1	0.0004	-1.0000	(0.0004, 4.0000)
2	0.0008	-1.0011	(0.0008, 4.0000)
3	0.0012	-1.0022	(0.0012, 4.0000)
6346	2.5384	-0.0064	(2.5384, 4.0000)
6347	2.5388	-0.0020	(2.5388, 4.0000)
6348	2.5388	0.0025	(2.5384, 2.5388)

搜索结束:
找到根：x = 2.5388, f(x) = 0.0025

在测试的过程中，发现最大迭代次数和步长和阈值需要相配合才出结果，步长大了可能就超出阈值永远也找不到根，迭代次数少了还没找到根就结束了，总之要找到一个误差较小的根，需要很多次运算，并且该方法的步长是固定的，也就意味着要想精度高，就得事先设置较小的步长，增加运算量

2.区间二分法

定理函数f(x)在[a,b]上单调连续，且f(a)f(b)<0，则方程f(x)=0在区间[a,b]上有且仅有一个实根x。
二分法的基本思想
将有根的区间二分为两个小区间，然后判断根在那个小区间，舍去无根的小区间，而把有根的小区间再一分为二，再判断根属于哪个更小的区间，如此反复，直到求出满足精度要求的近似根。
区间二分法的步骤：

选择初始区间： 首先，需要选择一个包含根的初始区间 [a, b]，其中 f(a) 和 f(b) 的符号不同，即 f(a) * f(b) < 0。这确保了根位于区间 [a, b] 内。
迭代过程： 利用以下迭代过程来逼近根的位置：

计算区间的中点： c = (a + b) / 2。
计算函数在中点的值： f(c)。
如果 f(c) 接近零（可以设置一个小的阈值来判断），则认为 c 是一个根，停止搜索。
否则，根据 f(a) 和 f(c) 的符号（或 f(b) 和 f(c) 的符号）来确定新的区间：
- 如果 f(a) * f(c) < 0，则根位于区间 [a, c] 内，更新 b = c。
- 如果 f(b) * f(c) < 0，则根位于区间 [c, b] 内，更新 a = c。

重复迭代： 重复步骤 2，直到找到满足条件的根或达到搜索的终止条件。通常，可以设置一个最大迭代次数或一个收敛阈值，以确保算法终止。
输出结果： 当算法终止时，得到的 c 就是近似的根，可以将其作为根的估计值输出。

二分过程结束原则：
先给定精度要求ε(绝对误差限)，
(1)事先由ε估计出二分的最小次数 k ，取 x≈xk

(2)当|bk+1 – ak+1|< ε时结束二分计算，取x≈xk；

区间二分法的优点是简单，收敛速度与比值为1/2的等比级数相同，它的局限性是只能用于求单根，不能求重根和复根。

测试

import math
def interval_bisection(f, a, b, tol=1e-6, max_iter=100):
    if f(a) * f(b) >= 0:
        raise ValueError("初始区间 [a, b] 不满足 f(a) * f(b) < 0 的条件")

    iteration = 0
    print("| 迭代次数 | 当前区间 [a, b] | f(a) * f(b) | f(b) * f(c) |")

    while iteration < max_iter:
        c = (a + b) / 2
        f_c = f(c)
        print(
            f"| {iteration + 1} | [{a}, {b}] | { '+' if f(a) * f(b)>0 else '-'} | {'+' if f(a) * f(c)>0 else '-'} |")

        if abs(f_c) < tol:
            return c  # 找到根，迭代结束

        if f(a) * f_c < 0:
            b = c
        else:
            a = c

        iteration += 1

    return (a + b) / 2  # 返回近似的根

def f(x):
    return x**3-4*x*math.log(x+2)-1

# 设置初始区间和收敛阈值
a = 0
b = 4
tolerance = 1e-6

# 调用区间二分法函数来寻找根
root = interval_bisection(f, a, b, tolerance)

print("找到的根的估计值为:", root)

迭代次数	当前区间 [a, b]	f(a) * f(b)	f(b) * f©
1	[0, 4]	-	+
2	[2.0, 4]	-	-
3	[2.0, 3.0]	-	+
4	[2.5, 3.0]	-	-
5	[2.5, 2.75]	-	-
6	[2.5, 2.625]	-	-
7	[2.5, 2.5625]	-	+
8	[2.53125, 2.5625]	-	-
9	[2.53125, 2.546875]	-	-
10	[2.53125, 2.5390625]	-	+
11	[2.53515625, 2.5390625]	-	+
12	[2.537109375, 2.5390625]	-	+
13	[2.5380859375, 2.5390625]	-	+
14	[2.53857421875, 2.5390625]	-	-
15	[2.53857421875, 2.538818359375]	-	-
16	[2.53857421875, 2.5386962890625]	-	-
17	[2.53857421875, 2.53863525390625]	-	-
18	[2.53857421875, 2.538604736328125]	-	-
19	[2.53857421875, 2.5385894775390625]	-	-
20	[2.53857421875, 2.5385818481445312]	-	-
21	[2.53857421875, 2.5385780334472656]	-	+
22	[2.538576126098633, 2.5385780334472656]	-	+
23	[2.538577079772949, 2.5385780334472656]	-	-

找到的根的估计值为: 2.5385775566101074

二、迭代法

基本迭代法或简单迭代法又称不动点迭代法、Picard迭代法等，为简单起见一般常称之为送代法。
介值定理 函数 f(x)在[a,b]上单调连续，且f(a)f(b)<0，则方程 f(x)=0 在区间[a,b]上有且仅有一个实根 x*。
**微分中值定理 **如果函数f(x)在[a,b]连续，在(a,b)可微，则在(a,b)内至少有一点ξ存在，使f(b)-f(a) =f’(ξ)(b-a) a<ξ<b

1.迭代原理

选代法是一种逐次逼近的方法，它是用某个固定公式反复校正根的近似值，使之逐步精确，最后得到满足精度要求的结果。迭代过程实质上是一个逐步显式化的过程数值计算方法第3版马东升_13878688_upscayl_4x_RealESRGAN_General_x4_v3.jpg

2.迭代的收敛性

定理（充分性条件）(全局收敛定理) 设函数 φ (x) 在[a, b]上连续，且
(1) 对 x∈[a, b],有 φ (x) ∈[a, b] (映内性)
(2) 存在0< L < 1，使对任意 x ∈[a, b]有|φ′(x)| ≤ L < 1 (压缩性)
则方程x= φ(x)在 [a, b]上的根x存在且唯一；
对初值x0 ∈[a, b] ,迭代过程 xk+1= φ (xk)均收敛于方程的根x。L越小收敛得越快。
一般不用全局收敛定理
**推论在定理的条件下，有误差估计式

定理设在区间[a, b]上方程 x= φ(x) 有根x,且对一切x∈[a, b]都有|φ′(x)| ≥ 1，则对于该区间上任意x0(≠x), 迭代公式xk+1= φ(xk)一定发散。
:::info
**局部收敛定理 **若在x = g(x)的根x*的某邻域内，有:
g’(x)连续，且存在0< L<1使得|g’(x)|≤L，（压缩性)
则对这邻域内的任意x，简单迭代法均收敛.
注意，这个定理要求找到根的邻域，因此先要使用零点存在性定理判别这个邻域
:::
:::info
迭代法的终点判断：只要相邻两次迭代值的偏差充分小，就能保证迭代值足够准确，因而用|xk- xk-1|控制迭代过程的结束。
:::

3.迭代的加速

加权法

埃特金迭代

斯蒂芬森迭代

迭代的计算步骤

确定 f(x)=0 的等价形式 x=φ(x) 选初值x0，判断收敛性|φ′(x0)|<1。
由公式x1 = φ(x0)计算x1
如果|x1-x0|≤ε则停止计算，取x* = x1 ；否则令 x0 = x1，重复步骤2 和 3，直到计算停止。

数值计算方法第3版马东升_13878688_upscayl_4x_RealESRGAN_General_x4_v3.jpg

测试

import math

#f(x)=x**3-4*x*math.log(x+2)-1
#φ′(x0)=0.43<1
def phi(x):
    return (4*x*math.log(x+2)+1)**(1/3)

def root_by_iteration(x0, epsilon):
    x1 = phi(x0)
    while abs(x1 - x0) > epsilon:
        x0 = x1
        x1 = phi(x0)
    return x1

# 设置初始值和精度
initial_value = 2.5
precision = 1e-10

# 调用函数进行迭代计算
root = root_by_iteration(initial_value, precision)

# 输出结果
print("迭代根:", root)

迭代根: 2.5385775512521187

加权法

import math

def phi(x):
    return (4*x*math.log(x+2)+1)**(1/3)

def weighted_iteration(x0, epsilon, max_iterations, weight):
    for _ in range(max_iterations):
        x1 = phi(x0)
        x0 = weight * x1 + (1 - weight) * x0
        if abs(x1 - x0) < epsilon:
            return x0
    return x0

# 设置初始值、精度、最大迭代次数和权重
initial_value = 2.5
precision = 1e-10
max_iterations = 100
weight = 0.5

# 调用函数进行迭代计算
root = weighted_iteration(initial_value, precision, max_iterations, weight)

# 输出结果
print("加权根:", root)

加权根: 2.5385775511274824

埃特金迭代

import math

def phi(x):
    return (4*x*math.log(x+2)+1)**(1/3)

def aitken_iteration(x0, epsilon, max_iterations):
    x1 = phi(x0)
    x2 = phi(x1)
    for _ in range(max_iterations):
        if abs(x2 - x1) < epsilon:
            return x2
        else:
            x0_temp = x0
            x1_temp = x1
            x2_temp = x2
            x0 = x2_temp - ((x2_temp - x1_temp)**2) / (x2_temp - 2*x1_temp + x0_temp)
            x1 = phi(x0)
            x2 = phi(x1)
    return x2

# 设置初始值、精度和最大迭代次数
initial_value = 2.5
precision = 1e-10
max_iterations = 100

# 调用函数进行迭代计算
root = aitken_iteration(initial_value, precision, max_iterations)

# 输出结果
print("埃特金迭代根:", root)

埃特金迭代根: 2.5385775513206923

斯蒂芬森迭代

import math

def phi(x):
    return (4*x*math.log(x+2)+1)**(1/3)

def steffensen_iteration(x0, epsilon, max_iterations):
    for _ in range(max_iterations):
        x1 = phi(x0)
        x2 = phi(x1)
        x0_temp = x0
        x0 = x0 - ((x1 - x0)**2) / (x2 - 2*x1 + x0)
        if abs(x0 - x0_temp) < epsilon:
            return x0
    return x0

# 设置初始值、精度和最大迭代次数
initial_value = 2.5
precision = 1e-10
max_iterations = 100

# 调用函数进行迭代计算
root = steffensen_iteration(initial_value, precision, max_iterations)

# 输出结果
print("斯蒂芬森迭代根:", root)

斯蒂芬森迭代根: 2.5385775513097064

三、牛顿迭代法(切线法)

牛顿迭代法的步骤：

给定一个初始近似值 x0。
计算函数在 x0 处的函数值 f(x0) 和导数值 f’(x0)。
使用牛顿迭代公式进行迭代：x1 = x0 - f(x0) / f’(x0)。
更新 x0 的值为 x1。
重复步骤 2 - 4，直到满足停止条件。

停止条件可以是以下之一：

达到预定的迭代次数。
函数值 f(x) 的绝对值小于或等于某个预定的容差值。
迭代过程中 x 的变化量小于或等于某个预定的容差值。

局部收敛定理：

若在f(x)=0的根x的某邻域内，有f"(x)连续，且f’(x)≠0,
则对这个邻域内的任意初值x0，牛顿迭代法均收敛。

全局收敛定理：

若f(x)= 0的根x∈[a,b]，且任意x ∈[a,b]，且f’(x)≠0，f"(x*）≠ 0,
则对[a,b]内，能使f(x0)*f"(x0)>0的初值x0，均能使牛顿迭代法收敛。

测试

import math

def f(x):
    return x**3 - 4*x*math.log(x+2) - 1

def f_prime(x):
    return 3*x**2 - 4*(math.log(x+2) + x/(x+2))

def newton_method(f, f_prime, x0, tolerance):
    x = x0
    while True:
        fx = f(x)
        if abs(fx) < tolerance:
            break
        x = x - fx / f_prime(x)

    return x

x0 = 2.5
tolerance = 1e-10

root = newton_method(f, f_prime, x0, tolerance)
print("牛顿迭代法的根:", root)

牛顿迭代法的根: 2.5385775513099365

四、弦截法(割线法)

牛顿迭代法的优点是收敛速度快，但它有个明显的缺点就是每送代一次都要计算f’(xk)，如果函数f(x)比较复杂，计算f’(xk)可能很麻烦，尤其当|f’(xk)|很小时，计算需很精确，否则误差会很大。为了避免计算导数，改用差商代替导数，这就是弦截法。弦截法又分单点弦法和双点弦法。双点弦法的收敛速度低于牛顿迭代法，但又高于不动点迭代法，因此是工程计算中常用的方法。

1.单点弦法

Snipaste_2023-09-16_15-36-14_upscayl_4x_RealESRGAN_General_x4_v3.png
这个公式的几何意义是过两个点作弦，这个弦和x轴的交点即是根的新的近似值，因为弦的一个端点(x0,f(x0))始终不变，只有另一个端点变动，所以这种方法称为单点弦法。

单点弦法的步骤

选择初始点：选择一个适当的初始点作为迭代的起点，记作 x0。
计算函数值：计算函数在初始点处的值，即 f(x0)。
选择第二个点：选择一个与初始点不同的第二个点作为迭代的终点，记作 x1。
计算函数值：计算函数在第二个点处的值，即 f(x1)。
计算斜率：计算通过点 (x0, f(x0)) 和 (x1, f(x1)) 的直线的斜率。斜率的计算公式为：slope = (f(x1) - f(x0)) / (x1 - x0)
计算下一个点：使用斜率来估计下一个点的位置。根据直线的方程，下一个点的 x 坐标可以通过以下公式计算：x2 = x1 - f(x1) / slope
更新点的位置：将 x1 更新为 x2，并将 x2 作为下一次迭代的终点。
判断终止条件：重复步骤 4-7，直到满足终止条件。终止条件可以是达到预定的迭代次数或函数值的精度达到要求。
输出结果：当满足终止条件时，输出最终的近似根。

测试

import math

def f(x):
    return x**3 - 4*x*math.log(x+2) - 1

def single_point_secant_method(f, x0, precision):
    x1 = x0 + 0.01  # 选择初始终点 x1，初始斜率不为零，避免除零错误
    while abs(f(x1)) > precision:
        slope = (f(x1) - f(x0)) / (x1 - x0)
        x2 = x1 - f(x1) / slope
        x0 = x1
        x1 = x2
    return x1


# 调用单点弦法求解方程
x0 = 2.5
precision = 1e-10
root = single_point_secant_method(f, x0, precision)

print("单点弦法的根:", root)

单点弦法的根: 2.5385775513097566

2.双点弦法

Snipaste_2023-09-16_15-43-18_upscayl_4x_RealESRGAN_General_x4_v3.png
双点弦法在计算xk+1时要用到前面两点的值xk和xk-1，这种迭代法称两步法，使用这类
方法，在计算前必须先提供两个开始值x0和x1。
双点弦法的几何意义如图所示。先过B0(x0,f(x0))点和B1(x1,f(x1))点作弦，这个弦（或其延长线）和x轴的交点x2，即是根的新的近似值，再过B1点和B2(x2,f(x2))点作弦求出x3，依此类推，当收敛时可求出满足精度要求的xk。

双点弦法的步骤

选择初始点：选择两个初始点x0和x1，使得f(x0)和f(x1)有不同的符号。这两个初始点将确定初始的割线。
计算割线与x轴的交点：根据割线的斜率，使用线性插值的方法计算割线与x轴的交点，得到新的近似根x2。
更新初始点：将x1的值赋给x0，并将x2的值赋给x1。
检查终止条件：如果满足终止条件（例如，达到了指定的精度或迭代次数），则停止迭代并返回近似根。
重复步骤2至4：使用更新后的初始点，重复步骤2至4，直到满足终止条件。

测试

import math

def f(x):
    return x**3 - 4*x*math.log(x+2) - 1

def secant_method(f, x0, x1, precision):
    while abs(f(x1)) > precision:
        x2 = x1 - f(x1) * (x1 - x0) / (f(x1) - f(x0))
        x0 = x1
        x1 = x2
    return x1

x0 = 2
x1 = 3
precision = 1e-10

root = secant_method(f, x0, x1, precision)

print("双点弦法的根:", root)