python 代码 实现之RK4

简化版的RK4代码

# 定义函数f，表示微分方程dy/dx = y - 2x/y
def f(x, y):
    return y - 2 * x / y

def rk4_step02(f, x, y, h):
    k1 = f(x, y)
    k2 = f(x + h / 2, y + k1*h / 2)
    k3 = f(x + h / 2, y + k2*h / 2)
    k4 = f(x + h, y + k3*h)
    return y + (k1 + 2 * k2 + 2 * k3 + k4)*h/ 6

# 使用四阶龙格-库塔法求解微分方程的主函数
def solve_rk4(f, x0, y0, h, x_end):
    x = x0
    y = y0
    results = [(x, y)]

    while x < x_end:
        y = rk4_step02(f, x, y, h)
        x += h
        results.append((x, y))

    return results

# 设置初始条件和参数
x0 = 0.0
y0 = 1.0
h = 0.25
x_end = 1.0

# 使用solve_rk4函数求解微分方程
results = solve_rk4(f, x0, y0, h, x_end)

# 打印结果
print("x\t\t\ty")
for x, y in results:
    print(f"{x:.1f}\t\t{y:.6f}")

显示计算过程

# 定义函数f，表示微分方程dy/dx = y - 2x/y
def f(x, y):
    return y - 2 * x / y

# 定义RK4方法的一个步骤，使用经典四阶龙格-库塔法计算下一个y值
def rk4_step01(f, x, y, h):
    print(f"rk4_step01: 计算开始，当前x={x}, y={y}, h={h}")
    k1 = h * f(x, y)
    k2 = h * f(x + h / 2, y + k1 / 2)
    k3 = h * f(x + h / 2, y + k2 / 2)
    k4 = h * f(x + h, y + k3)
    print(f"rk4_step01: k1={k1}, k2={k2}, k3={k3}, k4={k4}")
    result = y + (k1 + 2 * k2 + 2 * k3 + k4) / 6
    print(f"rk4_step01: 下一个y值为={result}")
    return result

# 定义RK4方法的另一个步骤，与rk4_step01类似，但有细微差别
def rk4_step02(f, x, y, h):
    print(f"rk4_step02: 计算开始，当前x={x}, y={y}, h={h}")
    k1 = f(x, y)
    k2 = f(x + h / 2, y + k1*h / 2)
    k3 = f(x + h / 2, y + k2*h / 2)
    k4 = f(x + h, y + k3*h)
    print(f"rk4_step02: k1={k1}, k2={k2}, k3={k3}, k4={k4}")
    result = y + (k1 + 2 * k2 + 2 * k3 + k4)*h / 6
    print(f"rk4_step02: 下一个y值为={result}")
    return result

# 使用四阶龙格-库塔法求解微分方程的主函数
def solve_rk4(f, x0, y0, h, x_end):
    x = x0
    y = y0
    results = [(x, y)]
    print("solve_rk4: 开始求解微分方程")

    while x < x_end:
        print(f"当前迭代：x={x}, y={y}")
        y = rk4_step02(f, x, y, h)
        x += h
        results.append((x, y))
        print(f"迭代完成：x={x}, y={y}")
        print(f"***************************************************************************************")

    print("solve_rk4: 求解完成")
    return results

# 设置初始条件和参数
x0 = 0.0
y0 = 1.0
h = 0.25
x_end = 1.0

# 使用solve_rk4函数求解微分方程
results = solve_rk4(f, x0, y0, h, x_end)

# 打印结果
print("x\t\t\ty")
for x, y in results:
    print(f"{x:.1f}\t\t{y:.6f}")