ODE Python 시스템을 해결하기위한 Runge-Kutta 4

Question

ODE 시스템을 해결하기 위해 Runge-Kutta 4 용 코드를 작성했습니다.
1D ODE에서는 잘 작동하지만 x'' + kx = 0해결하려고 할 때 벡터 함수를 정의하는 데 문제가 있습니다.

하자 u1 = x및 u2 = x' = u1'다음 시스템 외모가 좋아 :

u1' = u2
u2' = -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 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()
    

미리 감사드립니다.

Answer 1

올바른 길을 가고 있지만 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()

Answer 2

모델은 이것입니다 : 여기에 이미지 설명을 입력하십시오

Langtangen의 책 Programming for Computations-Python에서 발췌.