Runge-Kutta程序绘图异常排查：数组填充还是值调用错误？

我写了一个用Runge-Kutta方法求解常微分方程（ODE）并对比绘图的程序，理想情况下两个图应该完全重合，但Runge-Kutta函数输出的图不对。想问是数组填充错了，还是新值调用有问题？

原程序代码 #

#Import correct libraries and extensions
import numpy as np
from scipy.integrate import odeint
from matplotlib import pyplot as plt
import warnings
def fxn():
    warnings.warn("deprecated", DeprecationWarning)
with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    fxn()

#Define conditions and store values
h=0.02 #Given values for ODE
x0=1
y0=5
xpoints=[x0] #array for storing X values
ypoints=[y0] #Array for storing Y values

#State equations and define functions
def dy_dx(y,x):
    return x / (np.exp(x) - 1) #Provided Equation
def RungeKuttaFehlberg(x,y):
    return x / (np.exp(x) - 1)
#Calculates k1-k4, x and y solutions
def RKFAlg(x0,y0,h):
    k1 = RungeKuttaFehlberg(x0,y0)
    k2 = RungeKuttaFehlberg(x0+(h/2),y0+((h/2)*k1))
    k3 = RungeKuttaFehlberg(x0+(h/2),y0+((h/2)*k2))
    k4 = RungeKuttaFehlberg(x0+h,y0+(h*k3))
    y1 = y0+(h/6)*(k1+(2*k2)+(2*k3)+k4)
    x1 = x0+h
    x1 = round(x1,2)
    print("Point (x(n+1),y(n+1)) =",(x1,y1))
    return((x1,y1)) #Returns as ordered pair

#Define range for number of calculations
for i in range(2000):
    print(f"Y{i+1}".format(i)) #Solution value format
    x0,y0 = RKFAlg(x0,y0,h) #Calls RKF Function
    xpoints.append(x0) #Saves values into array
    ypoints.append(y0)
    y0 = 1
    ODEy1 = odeint(dy_dx,y0,xpoints)

#Runge-Kutta Graph
plt.plot(xpoints,ypoints,'b:',linewidth = 1) #command to plot lines using various colors and widths
plt.suptitle("RKF Graph")
plt.xlabel("x Points")
plt.ylabel("y Points")
plt.show()

#ODE graph
plt.plot(xpoints,ODEy1,'g-',linewidth=1)
plt.suptitle("ODE Graph")
plt.xlabel("x Points")
plt.ylabel("y Points")
plt.show()

#Function for plotting RKF and ODE graph
plt.plot(xpoints,ODEy1,'g-',linewidth=2,label="ODE")
plt.plot(xpoints,ypoints,'b:',linewidth=3,label="Runge-Kutta")
plt.suptitle("ODE and RKF Comparison")
plt.legend(bbox_to_anchor=(.8,1),loc=0,borderaxespad=0)
plt.xlabel("X")
plt.ylabel("Y")
plt.show()

#Function for plotting the difference graph
diff = [] #array to store difference
for i in range(len(xpoints)):
    diff.append(ypoints[i]-ODEy1[i])
plt.plot(xpoints,diff)
plt.suptitle("Difference")
plt.xlabel("x Points")
plt.ylabel("RKF and ODE diff.")
plt.show()

错误分析与修正点 #

• 核心错误：循环中强制重置y0为1：每次用RKFAlg算出下一个y值后，立刻执行y0 = 1，直接覆盖了正确的迭代初始值，导致后续所有RKF计算都从y=1开始，完全偏离了正确的积分路径。删掉这行代码即可解决RKF曲线的根本问题。
• odeint调用时机错误：把odeint放在循环内部会重复计算，而且每次用的是被重置后的y0=1，应该在循环结束后，用原始初始条件y0=5来计算整个xpoints对应的ODE解。
• 冗余函数定义：RungeKuttaFehlberg和dy_dx完全一致，直接复用dy_dx就行，不用重复写。
• x值四舍五入的累积误差：每次把x1四舍五入到两位小数，长期循环会让xpoints的步长偏离设定的0.02，建议去掉round操作，保留浮点精度。
• 维度不匹配问题：odeint返回的是二维数组，计算差值时需要转成一维，用ODEy1.flatten()处理。

修正后的代码 #

import numpy as np
from scipy.integrate import odeint
from matplotlib import pyplot as plt
import warnings

# 忽略警告（保留原逻辑）
def fxn():
    warnings.warn("deprecated", DeprecationWarning)
with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    fxn()

# 初始条件与参数
h = 0.02
x0 = 1
y0 = 5
xpoints = [x0]
ypoints = [y0]

# 定义微分方程（复用这一个函数即可）
def dy_dx(y, x):
    return x / (np.exp(x) - 1)

# RKF4算法实现
def RKFAlg(x_current, y_current, step):
    k1 = dy_dx(y_current, x_current)
    k2 = dy_dx(y_current + (step/2)*k1, x_current + step/2)
    k3 = dy_dx(y_current + (step/2)*k2, x_current + step/2)
    k4 = dy_dx(y_current + step*k3, x_current + step)
    y_next = y_current + (step/6)*(k1 + 2*k2 + 2*k3 + k4)
    x_next = x_current + step
    return x_next, y_next

# 迭代计算RKF解
for i in range(2000):
    x0, y0 = RKFAlg(x0, y0, h)
    xpoints.append(x0)
    ypoints.append(y0)

# 用odeint计算ODE的参考解（初始条件用y=5）
ODEy1 = odeint(dy_dx, 5, xpoints)
ODEy1_flat = ODEy1.flatten()  # 转成一维数组方便后续计算

# 绘图对比
plt.figure(figsize=(12, 8))

# RKF与ODE对比图
plt.subplot(2, 1, 1)
plt.plot(xpoints, ODEy1_flat, 'g-', linewidth=2, label="ODE")
plt.plot(xpoints, ypoints, 'b:', linewidth=3, label="Runge-Kutta")
plt.title("ODE和RKF求解结果对比")
plt.legend()
plt.xlabel("X")
plt.ylabel("Y")

# 差值图
plt.subplot(2, 1, 2)
diff = [yp - odeyp for yp, odeyp in zip(ypoints, ODEy1_flat)]
plt.plot(xpoints, diff)
plt.title("RKF与ODE解的差值")
plt.xlabel("X")
plt.ylabel("差值")

plt.tight_layout()
plt.show()

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