为何RK4方法在双摆仿真中能量异常波动而非缓慢衰减？

我清楚RK4方法不具备能量守恒特性，正常情况下能量会伴随小幅振荡缓慢衰减，但在我的双摆仿真中没有出现这种常规表现，具体情况如下：

我采用RK4方法实现双摆仿真，用Tkinter做动画，关闭窗口后绘制能量随时间变化的曲线，结果如下：

（能量单位：10千焦，时间为任意单位）

从图中可见，能量均值看似稳定，但与初始值偏差较大，且存在大幅振荡。我认为这个结果对于RK4方法而言很反常，特寻求解释。

代码结构说明 #

• __init__函数：初始化仿真变量与Tkinter窗口
• animate方法：包含RK4求解逻辑与动画更新代码
• calculate_energy方法：计算总能量并加入能量列表self.E
• on_closing与run函数：仿真结束后返回能量列表用于绘图

完整代码 #

import math
import tkinter as tk
import numpy as np
import matplotlib.pyplot as plt

class DoublePendulum:
    def __init__(self,  length1, length2, angle1, angle2, p1, p2, gravity,m1,m2,dt):
        
        global root
        root= tk.Tk()
        root.protocol("WM_DELETE_WINDOW", self.on_closing)
        # Create a TKinter window
        self.WINDOWSIZE=700

        self.canvas = tk.Canvas(root, width=self.WINDOWSIZE, height=self.WINDOWSIZE)
        self.canvas.pack()
        
        color="red"
        x=self.WINDOWSIZE/2
        y=self.WINDOWSIZE/2
        # Set initial position and velocity of pendulum
        self.origin = (x, y)
        self.angle1 = math.radians(angle1)
        self.angle2 = math.radians(angle2)
        self.p1 = math.radians(p1)  
        self.p2 = math.radians(p2)
        self.length1 = length1
        self.length2 = length2
        self.gravity = gravity
        self.m1 = m1
        self.m2 = m2
        self.dt=dt
        self.E=[]
        
        # Calculate initial position of pendulum
        self.position1 = (x + length1 * math.sin(self.angle1), y + length1 * math.cos(self.angle1))
        self.position2 = (self.position1[0] + length2 * math.sin(self.angle2), self.position1[1] + length2 * math.cos(self.angle2))
        
        # Draw the pendulum arms and bobs
        self.arm1 = self.canvas.create_line(self.origin[0], self.origin[1], 
                                            self.position1[0], self.position1[1], 
                                            width=3,fill=color)
        self.arm2 = self.canvas.create_line(self.position1[0], self.position1[1], 
                                            self.position2[0], self.position2[1], 
                                            width=3,fill=color)
        self.bob1 = self.canvas.create_oval(self.position1[0]+25, self.position1[1]+25, 
                                             self.position1[0]-25, self.position1[1]-25,fill=color)
        self.bob2 = self.canvas.create_oval(self.position2[0]+15, self.position2[1]+15, 
                                             self.position2[0]-15, self.position2[1]-15,fill=color)
        
        
    
    def animate(self):
        # Calculate the new positions of the pendulum arms and bobs
        # based on the current positions, velocities, and gravity
        current_state = [self.angle1, self.p1, self.angle2, self.p2]
        #F_A is just here to make the code simpler using those variables
        F_A = lambda y1: [(y1[1] * y1[3] * math.sin(y1[0] - y1[2])) / (self.length1 * self.length2 * (self.m1 + self.m2 * math.sin(y1[0] - y1[2])**2)),
                (1 / (2 * self.length1**2 * self.length2**2 * (self.m1 + self.m2 * math.sin(y1[0] - y1[2])**2)**2)) * ((y1[1]**2 * self.m2 * self.length2**2) - (2 * y1[1] * y1[3] * self.m2 * self.length1 * self.length2 * math.cos(y1[0] - y1[2])) + (y1[3]**2 * (self.m1 + self.m2) * self.length1**2)) * math.sin(2 * (y1[0] - y1[2]))]

        
        F = lambda y2: [(y2[1] * self.length2 - y2[3] * self.length1 * math.cos(y2[0] - y2[2])) / (self.length1**2 * self.length2 * (self.m1 + self.m2 * math.sin(y2[0] - y2[2])**2)),
                -((self.m1 + self.m2) * self.gravity * self.length1 * math.sin(y2[0])) - A_1 + A_2,
                (y2[3] * (self.m1 + self.m2) * self.length1 - y2[1] * self.m2 * self.length2 * math.cos(y2[0] - y2[2])) / (self.m2 * self.length1 * self.length2**2 * (self.m1 + self.m2 * math.sin(y2[0] - y2[2])**2)),
                -(self.m2 * self.gravity * self.length2 * math.sin(y2[2])) + A_1 - A_2]
        
        
        A_1,A_2=F_A(current_state)        
        Y1 = [self.dt * x for x in F(current_state)]
        A_1,A_2=F_A(np.array(current_state)+np.array([0.5 * f_i for f_i in Y1]))
        Y2=[self.dt * x for x in F(np.array(current_state)+np.array([0.5 * f_i for f_i in Y1]))]
        A_1,A_2=F_A(np.array(current_state)+np.array([0.5 * f_i for f_i in Y2]))
        Y3=[self.dt * x for x in F(np.array(current_state)+np.array([0.5 * f_i for f_i in Y2]))]
        A_1,A_2=F_A(np.array(current_state)+np.array(Y3))
        Y4=[self.dt * x for x in F(np.array(current_state)+(Y3))]
        Y=(np.array(Y1)+np.array([2 * x for x in Y2])+np.array([2 * x for x in Y3])+np.array(Y4))
        [self.angle1, self.p1, self.angle2, self.p2] = np.array(current_state) + np.array([(1/6) * x for x in Y])
        
        self.position1 = (self.origin[0] + self.length1 * math.sin(self.angle1), self.origin[1] + self.length1 * math.cos(self.angle1))
        self.position2 = (self.position1[0] + self.length2 * math.sin(self.angle2),self.position1[1] + self.length2 * math.cos(self.angle2))

        
        # Update the positions of the pendulum arms and bobs
        self.canvas.coords(self.arm1, self.origin[0], self.origin[1], self.position1[0], self.position1[1])
        self.canvas.coords(self.arm2, self.position1[0], self.position1[1], self.position2[0], self.position2[1])
        self.canvas.coords(self.bob1, self.position1[0]-25, self.position1[1]-25, self.position1[0]+25, self.position1[1]+25)
        self.canvas.coords(self.bob2, self.position2[0]-25, self.position2[1]-25, self.position2[0]+25, self.position2[1]+25)
        

        self.calculate_energy()
        # Call the animate function again after a short delay
        self.canvas.after(8, self.animate)
        
        

        
    def calculate_energy(self):
        #the calculus necessary to plot the energy in real time
        
        # Calculate the y positions 
        y1 = self.position1[1]
        y2 = self.position2[1]
        
        #calculate omegas, the angles derivative
        current_state = [self.angle1, self.p1, self.angle2, self.p2]
        Fh = lambda y: [(y[1] * self.length2 - y[3] * self.length1 * math.cos(y[0] - y[2])) / (self.length1**2 * self.length2 * (self.m1 + self.m2 * math.sin(y[0] - y[2])**2)),
                (y[3] * (self.m1 + self.m2) * self.length1 - y[1] * self.m2 * self.length2 * math.cos(y[0] - y[2])) / (self.m2 * self.length1 * self.length2**2 * (self.m1 + self.m2 * math.sin(y[0] - y[2])**2))]
        [omega1,omega2]=Fh(current_state)
        
        #calculate velocities of the masses
        vx1 = self.length1 * omega1 * np.cos(self.angle1)
        vy1 = self.length1 * omega1 * np.sin(self.angle1)
        vx2 = vx1 + self.length2 * omega2 * np.cos(self.angle2)
        vy2 = vy1 + self.length2 * omega2 * np.sin(self.angle2)
        
        # Calculate the kinetic energy of each mass
        K1 = 0.5 * self.m1 * (vx1**2 + vy1**2)
        K2 = 0.5 * self.m2 * (vx2**2 + vy2**2)
        
        # Calculate the potential energy of each mass
        U1 = self.m1 * self.gravity * y1
        U2 = self.m2 * self.gravity * y2
        
        # Add the kinetic and potential energy for each mass to get the total energy
        E = K1 + K2 + U1 + U2
        
        #add the calculated energy
        self.E.append(E/10000)

    
    def on_closing(self):
        #make the prog return self.E on closing
        root.return_value = self.E
        # Destroy the main window
        root.destroy()

    def run(self):
        # Start the main loop of TKinter
        self.animate()
        root.mainloop()
        return root.return_value

# Create a DoublePendulum object and start the animation
DoublePendulum=DoublePendulum(150, 150, 100, 250, 10, 15, 9.81, 10, 10, 0.1)
E=DoublePendulum.run()
# Define the time array
dt=0.1
t = [i*dt for i in range(len(E))]

# Plot the energy as a function of time
plt.plot(t, E)

# Set the plot title and axis labels
plt.title('Energy as a function of time')
plt.xlabel('Time')
plt.ylabel('Energy')

# Show the plot
plt.show()

已做检查 #

我核对了所用的双摆动力学方程，未发现错误。尝试提炼最小可复现示例但失败，因为问题核心在于整体仿真表现。所有测试环节均正常运行，仿佛这种能量异常表现是预期行为。

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