Runge-Kutta 4 用于求解 ODE 系统 Python答案

【问题标题】：Runge-Kutta 4 for solving systems of ODEs PythonRunge-Kutta 4 用于求解 ODE 系统 Python
【发布时间】：2020-08-26 22:19:21
【问题描述】：

我为 Runge-Kutta 4 编写了用于求解 ODE 系统的代码。
它适用于一维 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()

提前，谢谢。

【问题讨论】：

你用的是什么python版本？ 1/6 的评估结果是什么？
您能解释一下您对 RK4 步骤的组件化应用的动机吗？你知道广播在 numpy 数组的上下文中意味着什么吗？也就是说，当你添加一个向量和一个标量时会发生什么？
我使用的是 Python 3.8。 1/6因子来源于模型的推导。
您不能添加向量和标量。但我认为我不会那样做。函数 f 是一个数组，然后我为每个解添加 [j]。
但是你知道。当您选择一个组件时，您使u1 成为一个标量。在下一阶段，您将此标量添加到状态向量。将 RK4 步骤替换为 Euler 步骤，并考虑少量时间步骤的算法逻辑，定义状态向量的哪些组件，设置哪些组件，哪些结果有效，哪些由于输入不可用而无效。最后，您插入了一些不必要且错误的并发症。

标签： python differential-equations scientific-computing

【解决方案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()

【讨论】：

问题是'for j in range...' :) 谢谢。

【解决方案2】：

型号是这样的： enter image description here

来自 Langtangen 的书 Programming for Computations - Python。

【讨论】：