ode数值解+定义微分方程 dy/dx = -y + x + 1

import numpy as np
from scipy.optimize import fsolve


def f(x, y):
    """定义微分方程 dy/dx = -y + x + 1

    Args:
        x (float): 自变量
        y (float): 因变量

    Returns:
        float: 在点(x,y)处的导数值
    """
    return -y + x + 1


def exact_solution(x):
    """微分方程的解析解 y = e^(-x) + x

    Args:
        x (float/ndarray): 自变量值或数组

    Returns:
        float/ndarray: 对应x值的解析解
    """
    return np.exp(-x) + x


# 参数设置
h = 0.1  # 步长
x0 = 0.0  # 初始 x 值
y0 = 1.0  # 初始 y 值（需满足 y(0)=1）
n = 10  # 迭代次数
x_values = [x0 + i * h for i in range(n + 1)]  # 生成等间距x值数组


def euler_method():
    """欧拉法数值求解常微分方程

    Returns:
        list: 各节点处的y值列表

    Note:
        公式: y_{n+1} = y_n + h * f(x_n, y_n)
    """
    y = [y0]
    print("\n欧拉法:")
    for i in range(n):
        # 计算当前步的导数和下一个y值
        x_n, y_n = x_values[i], y[i]
        f_n = f(x_n, y_n)
        y_next = y_n + h * f_n
        y.append(y_next)
        print(f"迭代 {i}: x = {x_n:.1f}, y = {y_n:.8f}, f(x, y) = {f_n:.8f}, y_next = {y_next:.8f}")
    return y


def improved_euler_method():
    """改进欧拉法（Heun方法）数值求解常微分方程

    Returns:
        list: 各节点处的y值列表

    Note:
        分预测-校正两步：
        1. 预测步: y* = y_n + h*f(x_n, y_n)
        2. 校正步: y_{n+1} = y_n + h/2*[f(x_n,y_n) + f(x_{n+1},y*)]
    """
    y = [y0]
    print("\n改进欧拉法 (Heun):")
    for i in range(n):
        # 预测步
        x_n, y_n = x_values[i], y[i]
        f_n = f(x_n, y_n)
        y_star = y_n + h * f_n

        # 校正步
        f_star = f(x_values[i + 1], y_star)
        y_next = y_n + (h / 2) * (f_n + f_star)
        y.append(y_next)
        print(
            f"迭代 {i}: x = {x_n:.1f}, y = {y_n:.8f}, f(x, y) = {f_n:.8f}, y* = {y_star:.8f}, f(x+h, y*) = {f_star:.8f}, y_next = {y_next:.8f}")
    return y


def trapezoidal_method():
    """梯形法（隐式方法）数值求解常微分方程

    Returns:
        list: 各节点处的y值列表

    Note:
        需解隐式方程：
        y_{n+1} = y_n + h/2*[f(x_n,y_n) + f(x_{n+1},y_{n+1})]
        使用scipy.fsolve进行数值求解
    """
    y = [y0]
    print("\n梯形法:")
    for i in range(n):
        x_n, y_n = x_values[i], y[i]
        x_next = x_values[i + 1]

        # 定义隐式方程并求解
        func = lambda y_next: y_n + (h / 2) * (f(x_n, y_n) + f(x_next, y_next)) - y_next
        y_next = fsolve(func, y_n)[0]
        y.append(y_next)
        print(f"迭代 {i}: x = {x_n:.1f}, y = {y_n:.8f}, x_next = {x_next:.1f}, y_next = {y_next:.8f}")
    return y


def rk2_method():
    """二阶龙格-库塔法（中点法）数值求解常微分方程

    Returns:
        list: 各节点处的y值列表

    Note:
        计算步骤：
        1. k1 = f(x_n, y_n)
        2. k2 = f(x_n + h/2, y_n + h/2*k1)
        3. y_{n+1} = y_n + h*k2
    """
    y = [y0]
    print("\n二阶龙格-库塔法 (RK2):")
    for i in range(n):
        x_n, y_n = x_values[i], y[i]
        # 计算中间斜率
        f_n = f(x_n, y_n)
        x_mid = x_n + h / 2
        y_mid = y_n + (h / 2) * f_n
        f_mid = f(x_mid, y_mid)

        # 更新y值
        y_next = y_n + h * f_mid
        y.append(y_next)
        print(
            f"迭代 {i}: x = {x_n:.1f}, y = {y_n:.8f}, f(x, y) = {f_n:.8f}, x_mid = {x_mid:.2f}, y_mid = {y_mid:.8f}, f_mid = {f_mid:.8f}, y_next = {y_next:.8f}")
    return y


# 5. 三阶龙格-库塔法 (RK3)
def rk3_method():
    """
    使用三阶龙格-库塔法求解常微分方程的数值解

    全局参数说明:
    y0    - float: 初始条件y(x0)的值
    n     - int: 迭代总次数
    h     - float: 步长
    x_values - list[float]: 预先计算的离散x值序列
    f     - function: 微分方程函数，形式为f(x, y)

    返回值:
    list[float] - 各节点处的数值解列表，包含初始值y0
    """
    y = [y0]  # 初始化结果列表，存储初始值

    print("\n三阶龙格-库塔法 (RK3):")
    for i in range(n):
        x_n, y_n = x_values[i], y[i]

        # 计算三个斜率值
        k1 = f(x_n, y_n)  # 阶段1：当前点的斜率

        # 阶段2：中间点的斜率（步长中点处）
        k2 = f(x_n + h / 2, y_n + (h / 2) * k1)

        # 阶段3：调整后的中间点斜率（步长终点处）
        k3 = f(x_n + h, y_n - h * k1 + 2 * h * k2)

        # 使用加权平均计算下一时刻的y值
        y_next = y_n + (h / 6) * (k1 + 4 * k2 + k3)

        y.append(y_next)
        print(
            f"迭代 {i}: x = {x_n:.1f}, y = {y_n:.8f}, k1 = {k1:.8f}, k2 = {k2:.8f}, k3 = {k3:.8f}, y_next = {y_next:.8f}")

    return y


# 6. 四阶龙格-库塔法 (RK4)
def rk4_method():
    """
    使用四阶龙格-库塔法求解常微分方程初值问题
    依赖全局变量:
        y0: 初始y值
        n: 迭代次数
        x_values: 离散化的x值数组
        h: 步长
        f: 微分方程函数f(x,y)
    返回值:
        list: 各节点上的y值计算结果列表
    """
    y = [y0]
    print("\n四阶龙格-库塔法 (RK4):")
    for i in range(n):
        # 计算四个斜率系数
        x_n, y_n = x_values[i], y[i]
        k1 = f(x_n, y_n)
        k2 = f(x_n + h / 2, y_n + (h / 2) * k1)
        k3 = f(x_n + h / 2, y_n + (h / 2) * k2)
        k4 = f(x_n + h, y_n + h * k3)
        # 计算下一步的y值
        y_next = y_n + (h / 6) * (k1 + 2 * k2 + 2 * k3 + k4)
        y.append(y_next)
        print(
            f"迭代 {i}: x = {x_n:.1f}, y = {y_n:.8f}, k1 = {k1:.8f}, k2 = {k2:.8f}, k3 = {k3:.8f}, k4 = {k4:.8f}, y_next = {y_next:.8f}")
    return y


# 7. 隐式欧拉法
def implicit_euler_method():
    """
    使用隐式欧拉法求解常微分方程初值问题
    依赖全局变量:
        y0: 初始y值
        n: 迭代次数
        x_values: 离散化的x值数组
        h: 步长
        f: 微分方程函数f(x,y)
    返回值:
        list: 各节点上的y值计算结果列表
    """
    y = [y0]
    print("\n隐式欧拉法:")
    for i in range(n):
        # 构造隐式方程并求解
        x_n, y_n = x_values[i], y[i]
        x_next = x_values[i + 1]
        func = lambda y_next: y_n + h * f(x_next, y_next) - y_next
        y_next = fsolve(func, y_n)[0]
        y.append(y_next)
        print(f"迭代 {i}: x = {x_n:.1f}, y = {y_n:.8f}, x_next = {x_next:.1f}, y_next = {y_next:.8f}")
    return y


# 8. 显式亚当斯法 (二阶 Adams-Bashforth)
def adams_bashforth_2_method():
    """
    使用二阶Adams-Bashforth显式方法求解常微分方程初值问题
    依赖全局变量:
        y0: 初始y值
        x0: 初始x值
        n: 迭代次数
        x_values: 离散化的x值数组
        h: 步长
        f: 微分方程函数f(x,y)
    返回值:
        list: 各节点上的y值计算结果列表
    """
    y = [y0]
    # 使用显式欧拉法初始化第一步
    y1 = y0 + h * f(x0, y0)
    y.append(y1)
    print("\n显式亚当斯法 (二阶 Adams-Bashforth):")
    print(f"初始 (欧拉法): x = {x0:.1f}, y = {y0:.8f}, y1 = {y1:.8f}")
    # 使用二阶Adams公式迭代
    for i in range(1, n):
        x_n, y_n = x_values[i], y[i]
        f_n = f(x_n, y_n)
        f_prev = f(x_values[i - 1], y[i - 1])
        y_next = y_n + (h / 2) * (3 * f_n - f_prev)
        y.append(y_next)
        print(
            f"迭代 {i}: x = {x_n:.1f}, y = {y_n:.8f}, f(x, y) = {f_n:.8f}, f(x_prev, y_prev) = {f_prev:.8f}, y_next = {y_next:.8f}")
    return y


# 9. 显式亚当斯法 (四阶 Adams-Bashforth)
def adams_bashforth_4_method():
    """
    使用四阶Adams-Bashforth显式方法求解常微分方程初值问题
    依赖全局变量:
        y0: 初始y值
        n: 迭代次数
        x_values: 离散化的x值数组
        h: 步长
        f: 微分方程函数f(x,y)
    返回值:
        list: 各节点上的y值计算结果列表
    """
    y = [y0]
    # 使用RK4方法初始化前三个步骤
    temp_y = rk4_method()[:4]
    y.extend(temp_y[1:])
    print("\n显式亚当斯法 (四阶 Adams-Bashforth):")
    print(f"初始 (RK4): y = {temp_y[:3]}")
    # 应用四阶Adams公式计算后续步长
    x_n = x_values[2]
    y_n = y[2]
    f_n = f(x_n, y_n)
    f_n1 = f(x_values[1], y[1])
    f_n2 = f(x_values[0], y[0])
    y_next = y_n + (h / 24) * (55 * f_n - 59 * f_n1 + 37 * f_n2)
    y.append(y_next)
    print(
        f"迭代 2: x = {x_n:.1f}, y = {y_n:.8f}, f(x, y) = {f_n:.8f}, f(x-1, y-1) = {f_n1:.8f}, f(x-2, y-2) = {f_n2:.8f}, y_next = {y_next:.8f}")
    return y


# 10. 隐式亚当斯法 (二阶 Adams-Moulton)
def adams_moulton_2_method():
    """
    使用二阶Adams-Moulton隐式方法求解常微分方程初值问题
    依赖全局变量:
        y0: 初始y值
        x0: 初始x值
        n: 迭代次数
        x_values: 离散化的x值数组
        h: 步长
        f: 微分方程函数f(x,y)
    返回值:
        list: 各节点上的y值计算结果列表
    """
    y = [y0]
    # 使用显式欧拉法初始化第一步
    y1 = y0 + h * f(x0, y0)
    y.append(y1)
    print("\n隐式亚当斯法 (二阶 Adams-Moulton):")
    print(f"初始 (欧拉法): x = {x0:.1f}, y = {y0:.8f}, y1 = {y1:.8f}")
    # 构造并求解隐式方程
    for i in range(1, n):
        x_n, y_n = x_values[i], y[i]
        f_n = f(x_n, y_n)
        x_next = x_values[i + 1]
        func = lambda y_next: y_n + (h / 2) * (f(x_next, y_next) + f_n) - y_next
        y_next = fsolve(func, y_n)[0]
        y.append(y_next)
        print(
            f"迭代 {i}: x = {x_n:.1f}, y = {y_n:.8f}, f(x, y) = {f_n:.8f}, x_next = {x_next:.1f}, y_next = {y_next:.8f}")
    return y


# 11. 隐式亚当斯法 (四阶 Adams-Moulton)
def adams_moulton_4_method():
    """
    使用四阶Adams-Moulton隐式方法求解常微分方程

    步骤说明：
    1. 使用RK4方法计算前3个时间步的初始值
    2. 构建隐式方程并通过数值方法求解下一步值
    3. 采用单次迭代策略进行预测校正

    全局变量依赖说明：
    - y0: 初始条件值
    - x_values: 时间步序列
    - h: 步长
    - f: 微分方程右端函数

    返回值：
    list -- 各时间步的数值解列表，包含初始值和计算值
    """
    y = [y0]
    # 用 RK4 计算前三步
    # 获取前4个时间步的值(y0,y1,y2,y3)，但只保留后3个用于初始化
    temp_y = rk4_method()[:4]  # 得到 y0, y1, y2, y3
    y.extend(temp_y[1:])  # 添加 y1, y2, y3

    # 隐式方程求解核心部分
    print("\n隐式亚当斯法 (四阶 Adams-Moulton):")
    print(f"初始 (RK4): y = {temp_y[:3]}")
    # 准备前两步的函数值
    x_n = x_values[2]  # 当前时间点
    y_n = y[2]  # 当前状态量
    f_n = f(x_n, y_n)  # 当前导数值
    f_n1 = f(x_values[1], y[1])  # 前一步导数值

    # 构建隐式方程并求解
    x_next = x_values[3]
    # Adams-Moulton四阶公式：y_{n+1} = y_n + h/24*(9f(t_{n+1},y_{n+1}) + 19f_n -5f_{n-1})
    func = lambda y_next: y_n + (h / 24) * (9 * f(x_next, y_next) + 19 * f_n - 5 * f_n1) - y_next
    y_next = fsolve(func, y_n)[0]  # 使用数值方法求解非线性方程

    # 结果存储与输出
    y.append(y_next)
    print(
        f"迭代 2: x = {x_n:.1f}, y = {y_n:.8f}, f(x, y) = {f_n:.8f}, f(x-1, y-1) = {f_n1:.8f}, x_next = {x_next:.1f}, y_next = {y_next:.8f}")
    return y


# 执行所有方法并输出结果
# 数值方法字典：方法名称到对应计算函数的映射
# 包含常见的显式/隐式单步和多步方法，用于常微分方程数值解
methods = {
    "欧拉法": euler_method,
    "改进欧拉法": improved_euler_method,
    "梯形法": trapezoidal_method,
    "二阶 RK": rk2_method,
    "三阶 RK": rk3_method,
    "四阶 RK": rk4_method,
    "隐式欧拉法": implicit_euler_method,
    "显式 Adams (二阶)": adams_bashforth_2_method,
    "显式 Adams (四阶)": adams_bashforth_4_method,
    "隐式 Adams (二阶)": adams_moulton_2_method,
    "隐式 Adams (四阶)": adams_moulton_4_method
}

# 存储各方法计算结果（取x=0.5处的y值）
results = {}

# 遍历所有数值方法进行计算
# 每个方法返回完整y值序列，取末尾值作为x=0.3处的近似解
for name, method in methods.items():
    y_values = method()
    results[name] = y_values[-1]  # 取 y(0.5)

# 输出各方法在x=0.5处的数值解对比结果
# 格式化为表格形式，包含方法名称、近似值、解析解和绝对误差
print("\n结果对比 (x = 0.5):")
# 表头定义，左对齐并固定列宽
print(f"{'方法':<20} | {'y(0.5) 近似值':<15} | {'解析解':<15} | {'误差':<15}")
print("-" * 65)

# 计算精确解并与各数值方法的近似解进行对比
exact_y = exact_solution(0.5)
for name, y_approx in results.items():
    error = abs(y_approx - exact_y)
    # 格式化输出每行结果，保留8位小数
    print(f"{name:<20} | {y_approx:<15.8f} | {exact_y:<15.8f} | {error:<15.8f}")