Runge-Kutta 4 pour résoudre les systèmes d'ODE Python

Question

J'ai écrit du code pour Runge-Kutta 4 pour résoudre le système des ODE.
Cela fonctionne bien pour 1-D ODE, mais lorsque j'essaie de résoudre, x'' + kx = 0j'ai un problème en essayant de définir une fonction vectorielle:

Soit u1 = xet u2 = x' = u1', alors le système ressemble à:

u1' = u2
u2' = -k*u1

Si u = (u1,u2)et f(u, t) = (u2, -k*u1), alors nous devons résoudre:

u' = f(u, t)

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

Mon code entier est:

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

D'avance, merci.

Answer 1

Vous êtes sur la bonne voie, mais lorsque vous appliquez des méthodes d'intégration temporelle telles que RK à des ODE à valeurs vectorielles, on fait essentiellement exactement la même chose que dans le cas scalaire, juste avec des vecteurs.

Ainsi, vous sautez la for j in range(len(X_0))boucle et l'indexation associée et vous vous assurez de passer les valeurs initiales sous forme de vecteurs (tableaux numpy).

Également nettoyé tun peu l'indexation et stocké la solution dans une liste.

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

Le modèle est le suivant: entrez la description de l'image ici

Extrait du livre de Langtangen Programming for Computations - Python.