数值分析之数值微分

清华关陆版 p325

import numpy as np


# 两点公式（前向差分）
def two_point_forward(f, x, h):
    # 计算 f(x + h)
    f_x_plus_h = f(x + h)
    # 计算 f(x)
    f_x = f(x)
    # 计算前向差分结果
    result = (f_x_plus_h - f_x) / h

    # 打印前向差分计算细节
    print(f"前向差分计算细节:")
    print(f"f(x + h) = f({x} + {h}) = {f_x_plus_h}")
    print(f"f(x) = f({x}) = {f_x}")
    print(f"结果: (f(x + h) - f(x)) / h = ({f_x_plus_h} - {f_x}) / {h} = {result}")

    return result


# 两点公式（后向差分）
def two_point_backward(f, x, h):
    # 计算 f(x)
    f_x = f(x)
    # 计算 f(x - h)
    f_x_minus_h = f(x - h)
    # 计算后向差分结果
    result = (f_x - f_x_minus_h) / h

    # 打印后向差分计算细节
    print(f"后向差分计算细节:")
    print(f"f(x) = f({x}) = {f_x}")
    print(f"f(x - h) = f({x} - {h}) = {f_x_minus_h}")
    print(f"结果: (f(x) - f(x - h)) / h = ({f_x} - {f_x_minus_h}) / {h} = {result}")

    return result


# 两点中心差分
def two_point_central(f, x, h):
    # 计算 f(x + h)
    f_x_plus_h = f(x + h)
    # 计算 f(x - h)
    f_x_minus_h = f(x - h)
    # 计算两点中心差分结果
    result = (f_x_plus_h - f_x_minus_h) / (2 * h)

    # 打印两点中心差分计算细节
    print(f"两点中心差分计算细节:")
    print(f"f(x + h) = f({x} + {h}) = {f_x_plus_h}")
    print(f"f(x - h) = f({x} - {h}) = {f_x_minus_h}")
    print(f"结果: (f(x + h) - f(x - h)) / (2 * h) = ({f_x_plus_h} - {f_x_minus_h}) / (2 * {h}) = {result}")

    return result


# 三点公式（中心差分）
def three_point_central(f, x, h):
    # 计算 f(x + h)
    f_x_plus_h = f(x + h)
    # 计算 f(x - h)
    f_x_minus_h = f(x - h)
    # 计算三点中心差分结果
    result = (f_x_plus_h - f_x_minus_h) / (2 * h)

    # 打印三点中心差分计算细节
    print(f"三点中心差分计算细节:")
    print(f"f(x + h) = f({x} + {h}) = {f_x_plus_h}")
    print(f"f(x - h) = f({x} - {h}) = {f_x_minus_h}")
    print(f"结果: (f(x + h) - f(x - h)) / (2 * h) = ({f_x_plus_h} - {f_x_minus_h}) / (2 * {h}) = {result}")

    return result


# Richardson外推法改进三点中心差分公式
def richardson_extrapolation(f, x, h):
    # 使用步长 h 计算三点中心差分
    approx1 = three_point_central(f, x, h)
    # 使用步长 h/2 计算三点中心差分
    approx2 = three_point_central(f, x, h / 2)
    # 计算 Richardson 外推法结果
    result = (4 * approx2 - approx1) / 3

    # 打印 Richardson 外推法计算细节
    print(f"Richardson外推法计算细节:")
    print(f"三点中心差分 h = {h}: {approx1}")
    print(f"三点中心差分 h = {h / 2}: {approx2}")
    print(f"结果: (4 * approx2 - approx1) / 3 = (4 * {approx2} - {approx1}) / 3 = {result}")

    return result


# 五点公式（中心差分）
def five_point_central(f, x, h):
    # 计算 f(x - 2h)
    f_x_minus_2h = f(x - 2 * h)
    # 计算 f(x - h)
    f_x_minus_h = f(x - h)
    # 计算 f(x + h)
    f_x_plus_h = f(x + h)
    # 计算 f(x + 2h)
    f_x_plus_2h = f(x + 2 * h)
    # 计算五点中心差分结果
    result = (f_x_minus_2h - 8 * f_x_minus_h + 8 * f_x_plus_h - f_x_plus_2h) / (12 * h)

    # 打印五点中心差分计算细节
    print(f"五点中心差分计算细节:")
    print(f"f(x - 2h) = f({x} - 2 * {h}) = {f_x_minus_2h}")
    print(f"f(x - h) = f({x} - {h}) = {f_x_minus_h}")
    print(f"f(x + h) = f({x} + {h}) = {f_x_plus_h}")
    print(f"f(x + 2h) = f({x} + 2 * {h}) = {f_x_plus_2h}")
    print(
        f"结果: (f(x - 2h) - 8 * f(x - h) + 8 * f(x + h) - f(x + 2h)) / (12 * h) = ({f_x_minus_2h} - 8 * {f_x_minus_h} + 8 * {f_x_plus_h} - {f_x_plus_2h}) / (12 * {h}) = {result}")

    return result


# 示例函数
def f(x):
    return x * np.exp(x)


# 设置 x 的值和步长 h
x_value = 2.0
h_value = 0.1

# 设置 Richardson 外推法的步长
h_value_richardson = 0.2

# 打印不同差分方法的结果
print("\n两点公式（前向差分）结果:", two_point_forward(f, x_value, h_value))
print("\n\n")
print("\n两点公式（后向差分）结果:", two_point_backward(f, x_value, h_value))
print("\n\n")
print("\n两点公式(中心差分)结果:", two_point_central(f, x_value, h_value))
print("\n\n")
print("\n三点公式（中心差分）结果:", three_point_central(f, x_value, h_value))
print("\n\n")
print("\nRichardson外推法改进后的结果:", richardson_extrapolation(f, x_value, h_value_richardson))
print("\n\n")
print("\n五点公式（中心差分）结果:", five_point_central(f, x_value, h_value))