Runge-Kutta 4 để giải quyết các hệ thống ODEs Python

Question

Tôi đã viết mã cho Runge-Kutta 4 để giải quyết hệ thống ODE.
Nó hoạt động tốt cho ODE 1-D nhưng khi tôi cố gắng giải quyết, x'' + kx = 0tôi gặp sự cố khi cố gắng xác định một hàm vectơ:

Để u1 = xvà u2 = x' = u1', sau đó hệ thống trông giống như:

u1' = u2
u2' = -k*u1

Nếu u = (u1,u2)và f(u, t) = (u2, -k*u1), thì chúng ta cần giải quyết:

u' = f(u, t)

def f(u,t, omega=2):
    u, v = u
    return np.asarray([v, -omega**2*u])

Toàn bộ mã của tôi là:

import numpy as np

def 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, t

def 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()
    

Trước, cảm ơn bạn.

Answer 1

Bạn đang đi đúng hướng, nhưng khi áp dụng các phương pháp tích hợp thời gian chẳng hạn như RK cho các ODE có giá trị vectơ, về cơ bản người ta thực hiện điều tương tự như trong trường hợp vô hướng, chỉ với vectơ.

Do đó, bạn bỏ qua for j in range(len(X_0))vòng lặp và lập chỉ mục liên quan và bạn đảm bảo rằng bạn chuyển các giá trị ban đầu dưới dạng vectơ (mảng numpy).

Cũng làm sạch chỉ mục tmột chút và lưu trữ giải pháp trong một danh sách.

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()

Answer 2

Mô hình là thế này: nhập mô tả hình ảnh ở đây

Từ cuốn sách Lập trình cho tính toán - Python của Langtangen.