如何解决求解依赖于变量的 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 举报,一经查实,本站将立刻删除。