Python使用四阶龙格库塔法(RK4)求解常微分方程组代码校验咨询

原代码存在的错误 #

• 缩进错误：Python对缩进敏感，原代码中函数内部的变量定义、数组初始化等代码块没有正确缩进，运行会直接报语法错误
• 缺失依赖导入：没有导入numpy、matplotlib.pyplot两个依赖库
• 微分方程定义错误：和你需要求解的方程组不匹配：
  你需要求解的方程组为：

  x' = -x + sin(t)
  y' = -x - y
  对应状态向量x[0]为x(t)，x[1]为y(t)，导数向量第一个元素是x'(t)，第二个是y'(t)，原代码的导数计算完全不符合目标方程组

• 最后绘图缺少plt.show()调用，无法正常显示结果曲线

修正后可运行代码 #

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint # 用于结果验证

def GlucoseTT(x, t, params):
    a = params["a"]
    b = params["b"]
    c = params["c"]
    # 修正为匹配目标方程组的导数计算，此处a=b=c=1完全匹配你的需求
    xdot = np.array([-a*x[0] + np.sin(t), -b*x[0] - c*x[1]])
    return xdot

def RK4(f, x0, t0, tf, dt):
    t = np.arange(t0, tf, dt)
    nt = t.size
    nx = x0.size
    x = np.zeros((nx, nt))
    x[:,0] = x0
    for k in range(nt-1):
        k1 = dt * f(t[k], x[:,k])
        k2 = dt * f(t[k] + dt/2, x[:,k] + k1/2)
        k3 = dt * f(t[k] + dt/2, x[:,k] + k2/2)
        k4 = dt * f(t[k] + dt, x[:,k] + k3)
        dx = (k1 + 2*k2 + 2*k3 + k4)/6
        x[:,k+1] = x[:,k] + dx
    return x, t

# 定义问题参数
params = {"a":1, "b":1, "c":1}
f = lambda t, x : GlucoseTT(x, t, params)
x0 = np.array([0.75, 0])

# 求解ODE
t0 = 0
tf = 20 # 可自行改回100查看长周期变化，缩短区间更方便对比结果
dt = 0.1
x_rk4, t_rk4 = RK4(f, x0, t0, tf, dt)

# 用scipy官方ODE求解器验证结果
def ode_f(x, t):
    return GlucoseTT(x, t, params)
x_scipy = odeint(ode_f, x0, t_rk4)

# 绘图对比
plt.figure(figsize=(10,5))
plt.plot(t_rk4, x_rk4[0,:], "r", linewidth=3, label="RK4计算x(t)")
plt.plot(t_rk4, x_scipy[:,0], "k--", linewidth=1.5, label="scipy官方求解x(t)")
plt.plot(t_rk4, x_rk4[1,:], "b", linewidth=3, label="RK4计算y(t)")
plt.plot(t_rk4, x_scipy[:,1], "g--", linewidth=1.5, label="scipy官方求解y(t)")
plt.xlabel("t")
plt.ylabel("值")
plt.legend()
plt.grid(True)
plt.show()

验证结果说明 #

运行修正后代码可以看到，你实现的四阶龙格库塔法计算结果和scipy官方求解器的结果完全重合，说明RK4实现逻辑正确，求解结果符合要求。

内容的提问来源于stack exchange，提问作者Jojo98