目前，我正在嘗試解決 4 個耦合的 ODE，以穩定推車上的倒立擺。我用 Scipy 的 ODEINT 做這件事沒有問題，但是，我不能讓它與手動實作一起作業。這很可能是由于在代碼中的“模型”函式中進行了一些奇怪的資料格式化。
我嘗試了多種方法都無濟于事，因此我不會發布我的錯誤代碼，因為它們的范圍從在 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:\Users\Nikolai Lund Kühne\OneDrive - Aalborg Universitet\Uni\3. 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\3. 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 實作。
uj5u.com熱心網友回復：
我沒有完全除錯這個，你也可以將資料減少到預期發生的狀態。所以有些猜測。
Numpy 正以 Matlab 的風格發展。的構造格式K是一個行向量形狀的陣列，[[K1,K2,K3,K4]]。現在，任何形式的矩陣向量乘法，K@x，都有一個一維的結果。從數學上講，人們會期望標量或 1x1 矩陣[[u1]]。遵循 Matlab 哲學，它既不是，而是一個簡單的陣列u=[u1]。任何包含uinside 的進一步標量操作也將產生 1 元素陣列。將導數放在一起，這具有產生列向量的效果。現在對陣列的進一步操作有可能將其廣播到 4x4 矩陣形狀的陣列。4x4x4 形狀的張量是如何發生的我沒有跟進，但似乎很有可能。
轉載請註明出處，本文鏈接：https://www.uj5u.com/qiye/372135.html