算法报告1-解非线性方程(牛顿迭代法)

Python实现:
import numpy as np
import matplotlib.pyplot as plt

# 设置字体为 SimHei (黑体)
plt.rcParams['font.sans-serif'] = ['SimHei']

# 避免负号显示为方块
plt.rcParams['axes.unicode_minus'] = False

# 定义函数 f(x) 和 f'(x)
def f(x):
    return x ** 3 - np.sin(x) - 12 * x + 1


def df(x):
    return 3 * x ** 2 - np.cos(x) - 12


# 牛顿迭代法
def newton_iteration(x0, eps=1e-6, N=100):
    k = 0
    print("K----------Xk----------|Xk+1 - Xk|")
    print(f"{k}        {x0:.6f}")

    x_values = [x0]  # 用于可视化
    k += 1

    while k <= N:
        # 判断导数是否接近0,避免除零错误
        if abs(df(x0)) < 1e-6:
            print("迭代失败:导数为0。")
            return None, None

        # 计算x1
        x1 = x0 - f(x0) / df(x0)

        # 打印迭代信息
        print(f"{k}        {x1:.6f}       {abs(x1 - x0):.6f}")
        x_values.append(x1)  # 记录每次的迭代值

        # 判断是否收敛
        if abs(x1 - x0) < eps:
            print(f"找到解 x = {x1:.6f}, 经过 {k} 次迭代。")
            return x1, x_values

        # 更新x0
        x0 = x1
        k += 1

    # 达到最大迭代次数,迭代失败
    print("迭代失败:达到最大迭代次数N。")
    return None, None


# 可视化牛顿迭代法的过程
def visualize_newton(x_values, func):
    x_range = np.linspace(min(x_values) - 1, max(x_values) + 1, 400)
    y_values = func(x_range)

    plt.plot(x_range, y_values, label="f(x)")
    plt.axhline(0, color='black', linewidth=1)
    plt.axvline(0, color='black', linewidth=1)

    # 画出每次迭代的点
    for i, x in enumerate(x_values):
        plt.plot(x, func(x), 'ro')  # 画出当前点
        if i < len(x_values) - 1:
            plt.plot([x, x_values[i + 1]], [func(x), 0], 'g--')  # 画出垂直线
            plt.plot([x_values[i + 1], x_values[i + 1]], [0, func(x_values[i + 1])], 'g--') # 画出水平线

    plt.title("牛顿迭代法的迭代过程")
    plt.xlabel("x")
    plt.ylabel("f(x)")
    plt.grid(True)
    plt.legend()
    plt.show()


if __name__ == "__main__":
    # 输入初始值
    x0 = float(input("请输入初始值 x0: "))

    # 调用牛顿迭代法
    root, x_values = newton_iteration(x0)

    # 如果找到解,进行可视化
    if root is not None:
        visualize_newton(x_values, f)

函数 Python可视化图形代码:
import numpy as np
import matplotlib.pyplot as plt

# 定义函数 f(x) 和 f'(x)
def f(x):
    return x**3 - np.sin(x) - 12*x + 1

def f_prime(x):
    return 3*x**2 - np.cos(x) - 12

# 绘制 f(x)
x = np.linspace(-10, 10, 400)
y = f(x)

plt.plot(x, y, label='f(x)')
plt.axhline(0, color='black',linewidth=1)
plt.axvline(0, color='black',linewidth=1)
plt.legend()
plt.title('f(x) = x^3 - sin(x) - 12x + 1')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.grid(True)
plt.show()

 

posted @ 2024-09-24 22:06  Air_lwz  阅读(117)  评论(0)    收藏  举报