求解依赖于变量的 ODE，其解析表达式通过 scipy.integrate.solve_ivp 已知

问题描述

我正在尝试求解微分方程 dy/dt = f(t,y,k) 其中 k 是通过求解微分方程 dk/dx = g(t,k,x) (x 和 t是自变量）。下面是一个使用前向欧拉方法的示例。对于每个时间步，k 是通过使用其在前一个循环中的值求解其微分方程 get_dkdx 来计算的。

def forward_euler(f,t_eval,y0,k0):
    """ Forward Euler method
    :param f: function returning dy/dt 
    :param t_eval: time array
    :param y0: initial value of y
    :param k0: initial value of k """

    t = t_eval[0]
    y = y0
    k = k0
    x = np.linspace(0,1)
    tsol,ysol,ksol = [t],[y],[k]

    for i in range(len(t_eval) - 1):

        # Calculate the time step
        dt = t_eval[i+1] - t_eval[i]

        # Calculate the new value of y,t and k
        k = scint.odeint(lambda k_,x_: get_dkdx(t,k_,x_),k[0],x).T[0]  # solve with respect to x taking the
                                                                               # prevIoUs value of k as initial value
        y += dt * f(t,k)
        t += dt

        # Store the new values
        tsol.append(t)
        ysol.append(y)
        ksol.append(k)

    return np.array(ysol),np.array(ksol)

这是一个简单的例子：

def get_dydt(t,k):
    return np.sin(t) ** 2 * y - np.sum(k)  # random example


def get_dkdx(t,x):
    return k * x  # random example


import matplotlib.pyplot as plt

t1 = np.linspace(0,5,1001)
ysol1,ksol1 = forward_euler(get_dydt,t1,2,np.linspace(0,10))

figure,axes = plt.subplots(2)
axes[0].plot(t1,ysol1)
axes[1].plot(t1,ksol1)
axes[0].set_ylabel('y')
axes[1].set_ylabel('k')
axes[0].set_xlabel('t')
axes[1].set_xlabel('t')

但是，Forward Euler 方法不足以解决我的问题，我必须使用 scipy.integrate.solve_ivp 中的 Radau 或 RK45 方法。我该怎么做？

解决方法

暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！

如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@）