目前我正试图求解4个耦合的ODE，以稳定推车上的倒立摆。我用Scipy的ODEINT做这件事没有问题，但是，我不能用手动实现。这很可能是由于代码中的"model"函数中进行了一些奇怪的数据格式化。
我已经尝试了很多事情，但都无济于事，因此我不会发布我的错误代码，因为它们的大小在RK4方法中添加所有计算步骤时都不合适。
下面是我当前使用ODEINT的代码。我想问的是，是否有人能帮助我，这样函数"模型"就可以正确地制作出来，这样我就可以实现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:UsersNikolai Lund KühneOneDrive - Aalborg UniversitetUni3. semesterP3 - Dynamiske SystemermanualRK4.py", line 57, in <module>
a,b = RK4(model, 0, 12, x0, 1000)
File "C:UsersNikolai Lund KühneOneDrive - Aalborg UniversitetUni3. semesterP3 - Dynamiske SystemermanualRK4.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实现。
我没有完全调试它，您也可以将数据减少到预期发生的状态。所以有些猜测。
Numpy正以Matlab的风格出现。构造的K的格式是行向量[[K1,K2,K3,K4]]形状的阵列。现在，任何形式的矩阵向量乘法，K@x，都有一维的结果。在数学上，可以期望标量或1x1矩阵[[u1]]。根据Matlab的原理，它两者都不是，它是一个简单的数组u=[u1]。任何内部有u的进一步标量运算也将产生1元素数组。把导数放在一起，这就产生了一个列向量。现在，对阵列的进一步操作有可能将其广播到4x4矩阵形状的阵列。4x4x4形状的张量是如何发生的，我没有跟进，但这似乎很有可能。