欧拉法、改进欧拉法、二阶、三阶、四阶龙格-库塔法、隐式欧拉法、显式亚当斯法（二阶 Adams-Bashforth）和隐式亚当斯法（二阶 Adams-Moulton）。

大黄本 p281

import numpy as np
from scipy.optimize import fsolve

# 定义微分方程 f(x, y) = y - 2x/y
def f(x, y):
    return y - 2 * x / y

# 解析解 y = sqrt(2x + 1)
def exact_solution(x):
    return np.sqrt(2 * x + 1)

# 参数设置
h = 0.1  # 步长
x0 = 0.0  # 初始 x
y0 = 1.0  # 初始 y
n = 11    # 迭代次数
x_values = [x0 + i * h for i in range(n + 1)]  # x 的值: [0.0, 0.1, 0.2, 0.3]

# 1. 欧拉法
def euler_method():
    y = [y0]
    print("\n欧拉法:")
    for i in range(n):
        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

# 2. 改进欧拉法 (Heun 方法)
def improved_euler_method():
    y = [y0]
    print("\n改进欧拉法:")
    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

# 3. 二阶龙格-库塔法 (RK2)
def rk2_method():
    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_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

# 4. 三阶龙格-库塔法 (RK3)
def rk3_method():
    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)
        k2 = f(x_n + h / 2, y_n + (h / 2) * k1)
        k3 = f(x_n + h, y_n - h * k1 + 2 * h * k2)
        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

# 5. 四阶龙格-库塔法 (RK4)
def rk4_method():
    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_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

# 6. 隐式欧拉法
def implicit_euler_method():
    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]  # 使用 fsolve 求解
        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

# 7. 显式亚当斯法 (二阶 Adams-Bashforth)
def adams_bashforth_method():
    y = [y0]
    # 先用欧拉法计算 y1
    y1 = y0 + h * f(x0, y0)
    y.append(y1)
    print("\n显式亚当斯法 (二阶 Adams-Bashforth):")
    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)
        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

# 8. 隐式亚当斯法 (二阶 Adams-Moulton)
def adams_moulton_method():
    y = [y0]
    # 先用欧拉法计算 y1
    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

# 执行所有方法并输出结果
methods = {
    "欧拉法": euler_method,
    "改进欧拉法": improved_euler_method,
    "二阶 RK": rk2_method,
    "三阶 RK": rk3_method,
    "四阶 RK": rk4_method,
    "隐式欧拉法": implicit_euler_method,
    "显式 Adams": adams_bashforth_method,
    "隐式 Adams": adams_moulton_method
}

results = {}
for name, method in methods.items():
    y_values = method()
    results[name] = y_values[-1]  # 取 y(0.3)

# 输出结果对比
print("\n结果对比 (x = 0.5):")
print(f"{'方法':<15} | {'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)
    print(f"{name:<15} | {y_approx:<15.8f} | {exact_y:<15.8f} | {error:<15.8f}")