我為 Runge-Kutta 4 編寫了代碼來求解 ODE 系統。它對于一維 ODE 工作得很好，但是當我嘗試解決時，x'' + kx = 0我在嘗試定義向量函數時遇到了問題：讓u1 = x和u2 = x' = u1'，則系統如下所示：u1' = u2u2' = -k*u1如果u = (u1,u2)和f(u, t) = (u2, -k*u1)，那么我們需要解決：u' = f(u, t)def f(u,t, omega=2): u, v = u return np.asarray([v, -omega**2*u])我的整個代碼是：import numpy as npdef ode_RK4(f, X_0, dt, T): N_t = int(round(T/dt)) # Create an array for the functions ui u = np.zeros((len(X_0),N_t+1)) # Array u[j,:] corresponds to the j-solution t = np.linspace(0, N_t*dt, N_t + 1) # Initial conditions for j in range(len(X_0)): u[j,0] = X_0[j] # RK4 for j in range(len(X_0)): for n in range(N_t): u1 = f(u[j,n] + 0.5*dt* f(u[j,n], t[n])[j], t[n] + 0.5*dt)[j] u2 = f(u[j,n] + 0.5*dt*u1, t[n] + 0.5*dt)[j] u3 = f(u[j,n] + dt*u2, t[n] + dt)[j] u[j, n+1] = u[j,n] + (1/6)*dt*( f(u[j,n], t[n])[j] + 2*u1 + 2*u2 + u3) return u, tdef demo_exp(): import matplotlib.pyplot as plt def f(u,t): return np.asarray([u]) u, t = ode_RK4(f, [1] , 0.1, 1.5) plt.plot(t, u[0,:],"b*", t, np.exp(t), "r-") plt.show() def demo_osci(): import matplotlib.pyplot as plt def f(u,t, omega=2): # u, v = u Here I've got a problem return np.asarray([v, -omega**2*u]) u, t = ode_RK4(f, [2,0], 0.1, 2) for i in [1]: plt.plot(t, u[i,:], "b*") plt.show()
1 回答
ABOUTYOU
TA貢獻1812條經驗 獲得超5個贊
您走在正確的道路上，但是當將 RK 等時間積分方法應用于向量值 ODE 時，本質上與標量情況下執行的操作完全相同，只是使用向量。
因此，您可以跳過for j in range(len(X_0))循環和關聯的索引，并確保將初始值作為向量（numpy 數組）傳遞。
還清理了一點索引t并將解決方案存儲在列表中。
import numpy as np
def ode_RK4(f, X_0, dt, T):
N_t = int(round(T/dt))
# Initial conditions
usol = [X_0]
u = np.copy(X_0)
tt = np.linspace(0, N_t*dt, N_t + 1)
# RK4
for t in tt[:-1]:
u1 = f(u + 0.5*dt* f(u, t), t + 0.5*dt)
u2 = f(u + 0.5*dt*u1, t + 0.5*dt)
u3 = f(u + dt*u2, t + dt)
u = u + (1/6)*dt*( f(u, t) + 2*u1 + 2*u2 + u3)
usol.append(u)
return usol, tt
def demo_exp():
import matplotlib.pyplot as plt
def f(u,t):
return np.asarray([u])
u, t = ode_RK4(f, np.array([1]) , 0.1, 1.5)
plt.plot(t, u, "b*", t, np.exp(t), "r-")
plt.show()
def demo_osci():
import matplotlib.pyplot as plt
def f(u,t, omega=2):
u, v = u
return np.asarray([v, -omega**2*u])
u, t = ode_RK4(f, np.array([2,0]), 0.1, 2)
u1 = [a[0] for a in u]
for i in [1]:
plt.plot(t, u1, "b*")
plt.show()
添加回答
舉報
0/150
提交
取消