手动RK4解决IVP（数据格式化问题）的方法

Question

目前我正在尝试求解 4 个耦合 ODE 以稳定手推车上的倒立摆。我使用 Scipy 中的 ODEINT 没有问题，但是，我无法使其与手动实现一起使用。这很可能是由于在代码的 'model' 函数中进行了一些奇怪的数据格式化。

我尝试了多种方法都无济于事，因此我不会 post 我的错误代码，因为它们的范围是在 RK4 方法中添加所有计算步骤时大小不合适。

我当前使用 ODEINT 的代码在下面。我要问的是是否有人可以帮助我，以便正确制作函数 'model'，以便我可以实现 RK4 求解器（我可以毫无问题地为其他 ODE 求解）。

import numpy as np
from scipy.integrate import solve_ivp
from scipy import signal
g = 9.82
l = 0.281 
mc = 6.28
alpha = 0.4
mp = 0.175
t_start = 0.
t_end = 12.
tol = 10**(-1)
# Define A and B and the poles we want
A = np.array([[0., 1., 0., 0.], [(mc+mp)*g/(l*mc), 0., 0., (-alpha)/(l*mc)], [0., 0., 0., 1.], [(g*mp)/mc, 0., 0., (-alpha)/mc]])
B = np.array([[0.], [1./(l*mc)], [0.], [1./mc]])
Poles = np.array([complex(-1.,2.), complex(-1.,-2.), complex(-2.,1.), complex(-2.,-1.)])

# Determine K
signal = signal.place_poles(A, B, Poles)
K = signal.gain_matrix
# print(signal.computed_poles) # To verify if the computes poles are correct

# Define the model
def model(t,x):
    x1, x2, x3, x4 = x
    u = -np.matmul(K,x)
    dx1dt = x2
    dx2dt = (np.cos(x1.astype(float))*(u-alpha*x4-mp*l*x2**2*np.sin(x1.astype(float)))+(mc+mp)*g*np.sin(x1.astype(float)))/(l*(mc+mp*(1-np.cos(x1.astype(float))**2)))
    dx3dt = x4
    dx4dt = (u-alpha*x4-mp*l*x2**2*np.sin(x1.astype(float))+mp*g*np.sin(x1.astype(float))*np.cos(x1.astype(float)))/(mc+mp*(1-np.cos(x1.astype(float))**2))
    return np.array([dx1dt, dx2dt, dx3dt, dx4dt])

# Solve the system
N = 10000 # Number of steps
t = np.linspace(t_start, t_end, N)
t_span = (t_start, t_end)
x0 = np.array([0.2, 0., 0., 0.])
sol = solve_ivp(model,t_span,x0, t_eval=t, method='RK45')
index = np.argmin(sol.y[2,:]) # Max displacement from the origin
print(f' The biggest deviation from the origin is: {abs(sol.y[2, index])} meters.')

#This doesn't work
def RK4(fcn,a ,b ,y0 ,N):
    h = (b-a)/N
    x = a + np.arange(N+1)*h
    y = np.zeros((x.size,y0.size))
    y[0,:] = y0
    for k in range(N):
        k1 = fcn(x[k], y[k,:])
        k2 = fcn(x[k] + h/2, y[k,:] + h*k1/2)
        k3 = fcn(x[k] + h/2, y[k,:] + h*k2/2)
        k4 = fcn(x[k] + h, y[k,:] + h*k3)
        y[k+1,:] = y[k,:] + h/6*(k1 + 2*(k2 + k3) + k4)
    
    return x,y

a,b = RK4(model, 0, 12, x0, 1000)

这会产生以下错误：

runcell(0, 'C:/Users/Nikolai Lund Kühne/OneDrive - Aalborg Universitet/Uni/3. semester/P3 - Dynamiske Systemer/manualRK4.py')
 The biggest deviation from the origin is: 0.48256054833140316 meters.
Traceback (most recent call last):

  File "C:\Users\Nikolai Lund Kühne\OneDrive - Aalborg Universitet\Uni. semester\P3 - Dynamiske Systemer\manualRK4.py", line 57, in <module>
    a,b = RK4(model, 0, 12, x0, 1000)

  File "C:\Users\Nikolai Lund Kühne\OneDrive - Aalborg Universitet\Uni. semester\P3 - Dynamiske Systemer\manualRK4.py", line 53, in RK4
    y[k+1,:] = y[k,:] + h/6*(k1 + 2*(k2 + k3) + k4)

ValueError: could not broadcast input array from shape (4,4,4) into shape (4)

编辑 2：尝试手动实施 RK4 会导致一些奇怪的错误。

编辑 1：根据评论，代码现在使用 solve_ivp 实现。

Answer 1

我没有完全调试这个，你也可以将数据减少到预期发生的状态。所以一些猜测。

Numpy 正在模仿 Matlab 的风格。 K的构造格式是一个行向量形状的数组，[[K1,K2,K3,K4]]。现在任何形式的矩阵向量乘法 K@x 都有一个一维结果。从数学上讲，人们会期望标量或 1x1 矩阵 [[u1]]。遵循 Matlab 哲学，两者都不是，它是一个简单的数组 u=[u1]。任何包含 u 的进一步标量运算也将产生 1 元素数组。将导数放在一起，这具有产生列向量的效果。现在，对数组的进一步操作有可能将其广播到 4x4 矩阵形状的数组。 4x4x4 形张量是如何产生的我没有跟进，但似乎很有可能。

手动RK4解决IVP（数据格式化问题）的方法

Manual RK4 method for solving IVP (data formatting problem)

python

numpy

scipy

differential-equations

runge-kutta