y'=-50y+50x^2+2x,y(0)=1/3,0<=x<=1.步长h=0.001,迭代3次。分别用RK4法,Runge-Kutta-Fehlberg法求解。要求显示详细的演绎计算过程

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import ScalarFormatter


# --- 1. 问题定义 ---
# 定义常微分方程 y' = f(x, y)
def f(x, y):
    """ODE function: y' = -50y + 50x^2 + 2x"""
    return -50 * y + 50 * x ** 2 + 2 * x


# 定义精确解
def y_exact(x):
    """Exact solution: y(x) = (1/3)e^(-50x) + x^2"""
    return (1 / 3) * np.exp(-50 * x) + x ** 2


# --- 2. 数值求解器实现 ---
def rk4_solver(f, x0, y0, h, x_end):
    """四阶龙格-库塔法求解器"""
    x_points = [x0]
    y_points = [y0]

    xi, yi = x0, y0

    while xi < x_end:
        k1 = f(xi, yi)
        k2 = f(xi + 0.5 * h, yi + 0.5 * h * k1)
        k3 = f(xi + 0.5 * h, yi + 0.5 * h * k2)
        k4 = f(xi + h, yi + h * k3)

        yi = yi + (h / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4)
        xi = xi + h

        x_points.append(xi)
        y_points.append(yi)

    return np.array(x_points), np.array(y_points)


def rkf45_solver(f, x0, y0, h, x_end):
    """龙格-库塔-费尔伯格(RKF45)法求解器 (使用四阶公式)"""
    x_points = [x0]
    y_points = [y0]

    xi, yi = x0, y0

    # RKF45 系数
    c = [0, 1 / 4, 3 / 8, 12 / 13, 1, 1 / 2]
    a = [
        [],
        [1 / 4],
        [3 / 32, 9 / 32],
        [1932 / 2197, -7200 / 2197, 7296 / 2197],
        [439 / 216, -8, 3680 / 513, -845 / 4104],
        [-8 / 27, 2, -3544 / 2565, 1859 / 4104, -11 / 40]
    ]
    # 四阶解的系数
    b4 = [25 / 216, 0, 1408 / 2565, 2197 / 4104, -1 / 5, 0]

    while xi < x_end:
        k = np.zeros(6)
        k[0] = f(xi, yi)
        for i in range(1, 6):
            y_temp = yi + h * sum(a[i][j] * k[j] for j in range(i))
            k[i] = f(xi + c[i] * h, y_temp)

        # 使用四阶公式更新 y
        yi = yi + h * sum(b4[i] * k[i] for i in range(6))
        xi = xi + h

        x_points.append(xi)
        y_points.append(yi)

    return np.array(x_points), np.array(y_points)


# --- 3. 计算与可视化 ---
# 初始条件和参数
x0, y0 = 0.0, 1 / 3
h = 0.0001
x_end = 1.0

# 运行求解器
x_rk4, y_rk4 = rk4_solver(f, x0, y0, h, x_end)
x_rkf, y_rkf = rkf45_solver(f, x0, y0, h, x_end)

# 计算用于绘图的精确解
x_true = np.linspace(0, x_end, 500)
y_true = y_exact(x_true)

# 计算每个数值解时间点的误差
y_true_at_steps_rk4 = y_exact(x_rk4)
error_rk4 = np.abs(y_rk4 - y_true_at_steps_rk4)

y_true_at_steps_rkf = y_exact(x_rkf)
error_rkf = np.abs(y_rkf - y_true_at_steps_rkf)

# --- 4. 绘图 ---
# 设置中文字体以正确显示图例
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 10), sharex=True)

# 图1: 解的对比
ax1.plot(x_true, y_true, 'k-', label='精确解', linewidth=2.5)
ax1.plot(x_rk4, y_rk4, 'r--', label='RK4 解', linewidth=1.5)
ax1.plot(x_rkf, y_rkf, 'b:', label='RKF45 解', linewidth=1.5)
ax1.set_title('数值解与精确解对比', fontsize=16)
ax1.set_ylabel('y(x)', fontsize=12)
ax1.legend(fontsize=12)
ax1.grid(True)

# 放大显示初始阶段的细节
ax_zoom = ax1.inset_axes([0.1, 0.1, 0.4, 0.4])
ax_zoom.plot(x_true, y_true, 'k-')
ax_zoom.plot(x_rk4, y_rk4, 'r--')
ax_zoom.plot(x_rkf, y_rkf, 'b:')
ax_zoom.set_xlim(0, 0.1)
ax_zoom.set_ylim(0, 0.4)
ax_zoom.grid(True)
ax1.indicate_inset_zoom(ax_zoom, edgecolor="black")

# 图2: 误差对比
ax2.plot(x_rk4, error_rk4, 'r-', label='RK4 绝对误差')
ax2.plot(x_rkf, error_rkf, 'b-', label='RKF45 绝对误差')
ax2.set_yscale('log')  # 使用对数坐标轴，因为误差很小
ax2.set_title('绝对误差对比 (对数坐标)', fontsize=16)
ax2.set_xlabel('x', fontsize=12)
ax2.set_ylabel('绝对误差', fontsize=12)
ax2.legend(fontsize=12)
ax2.grid(True, which="both", ls="--")

# 格式化y轴刻度为科学计数法
formatter = ScalarFormatter(useMathText=True)
formatter.set_scientific(True)
formatter.set_powerlimits((-1, 1))
ax2.yaxis.set_major_formatter(formatter)

plt.tight_layout()
plt.show()