用RungeKutta法求解运动微分方程

Question

在我们的物理课上，我们必须建立一个阻尼扭转摆的模型

i扭摆的说明：

我们得出了这个运动方程：

其中，θ是角度，A是扭转参数，B是牛顿参数，C是斯托克斯参数，D是摩擦参数。我们还使用符号函数sgn，该函数根据当前与参考点的角度确定摆锤上作用力的方向

问题是，我无法使用Python中的龙格库塔方法解决它

我在MATLAB中使用Euler方法得到了一个有效的解决方案，它有一些缺陷，但它是有意义的

MATLAB代码：

function [theta, dtheta, epsilon] = drt(t, theta0, dtheta0, A, B, C, D)
    epsilon = zeros(1, length(t));
    theta = epsilon;
    dtheta = epsilon;

    theta(1) = theta0;
    dtheta(1) = dtheta0;


    epsilon(1) = A * alpha0 - B * dtheta^2 - C * dtheta - D;

    dt = t(2) - t(1);

    for i = 1 : (length(t) - 1)
       epsilon(i + 1)= - A * theta(i) - B * dtheta(i)^2 * sign(dtheta(i)) - C * dtheta(i) - D * sign(dtheta(i));
       dtheta(i + 1)= dtheta(i) + dt * epsilon(i); 
       theta(i + 1) = theta(i) + dt * dtheta(i);
    
    end
end

我们将此函数称为MATLAB函数，例如：

t = linspace(0, 10, 100);
theta0 = 90;
dtheta0 = 0;
A = 1;
B = 0.1;
C = 0.1;
D = 0.1;
[theta, dtheta, epsilon] = drt(t, theta0, dtheta0, A, B, C, D);

然后，我们可以在一张图表中绘制出θ和其他值，它向我们展示了扭转摆是如何被作用在其上的外力阻尼的

Python代码：

import numpy as np
import matplotlib.pyplot as plt


# Damping torsional pendulum
def drp(drp_alpha, drp_d_alpha, drp_params):
    a = drp_params["tors"]
    b = drp_params["newt"]
    c = drp_params["stok"]
    d = drp_params["fric"]

    result = a * drp_alpha - b * np.power(drp_d_alpha, 2) * np.sign(drp_d_alpha) - c * drp_d_alpha - d * np.sign(drp_d_alpha)

    return result


# Runge-Kutta 4th Order
# f - function DamRotPen
# x0 - initial condition
# t0 - initial time
# tmax - maximum time
# dt - sample time
def RG4(rg4_f, rg4_x0, rg4_t0, rg4_tmax, rg4_dt):
    # Time vector
    rg4_t = np.arange(rg4_t0, rg4_tmax, rg4_dt)
    # Time vector size
    rg4_t_sz = rg4_t.size

    # Initialize the array
    rg4_alpha = np.zeros(rg4_t_sz)
    # Initial value of the system
    rg4_alpha[0] = rg4_x0

    for k in range(rg4_t_sz - 1):
        k1 = dt * f(rg4_t[k], rg4_alpha[k])
        k2 = dt * rg4_f(rg4_t[k] + dt / 2, rg4_alpha[k] + k1 / 2)
        k3 = dt * rg4_f(rg4_t[k] + dt / 2, rg4_alpha[k] + k2 / 2)
        k4 = dt * rg4_f(rg4_t[k] + dt, rg4_alpha[k] + k3)

        rg4_d_alpha = (k1 + 2 * k2 + 2 * k3 + k4) / 6

        rg4_alpha[k + 1] = rg4_alpha[k] + rg4_d_alpha

    return rg4_alpha, rg4_t


# Parameters of the forces acting on the system
# tors - torsion parameter
# newt - Newton's parameter
# stok - Stokes' parameter
# fric - friction parameter
params = {"tors": 1, "newt": 0.1, "stok": 0.1, "fric": 0.1}

# Start parameters
alpha = 90
d_alpha = 0

# Initial time
t0 = 0
# Maximum time
tmax = 120
# Sample time
dt = 0.01

# Define DamRotPen function as 'f' using lambda
f = lambda t, alpha : drp(alpha, d_alpha, params)

# Try to solve this shit
alpha, t = RG4(f, alpha, t0, tmax, dt)

# Plot this shit
plt.plot(t, alpha, "r", "r", label="Position I guess")
plt.xlabel("Time t / s")
plt.grid()
plt.show()

绘制完这些值后，我们可以在图上看到，θ天空在一定时间后会发射火箭。我不知道我做错了什么，我几乎什么都试过了，所以我要你帮忙。（虽然我认为问题的一部分可能是我对如何实现龙格-库塔方法的误解，但可能是我的数学算错了，等等。）