问题描述
嘿,我有一个模型可以提供不同的返回值。但也是一个不是微分 (z) 的值
def model(t,f):
x = f[0]
y = f[1]
dx_dt = 2x*y
dy_dt = x**3
z = 2*x
return dx_dt,dy_dt,z
我的求解器可以求解这个方程,并在各自的时间值处给出 x 和 y。
t = np.linspace(0,10,100)
f0 = [2,1,0]
result = solve_ivp(model,[np.min(t),np.max(t)],f0,t_eval=t)
但现在我还想要我的 z 解决方案,该解决方案不应由求解器集成。有没有可能只对求解器的前 2 个值进行积分?但不是第三个?
解决方法
你有 ODE 函数可以写成这样的情况
def model(t,u):
v = f(t,u)
du = g(t,u,v) # v may contain components of du that get simply copied
return du
并且您还对数量 z=h(t,v)
感兴趣。所有变量也可以是元组或向量。目前,函数代表代码块,但在许多情况下,它们也可以很容易地作为函数分离出来。因此,第一个变体就是这样做,以便 ODE 函数具有最少的额外功能,并且值 z
在最后根据解值计算。
另一种变体,如果将代码块转换为单独的函数似乎不是最优的,您也可以按照问题构造模型函数,
def model(t,v) # v may contain components of du that get simply copied
z = h(t,v)
return du,z
然后在求解器的调用中使用 lambda 表达式来分离导数向量
result = solve_ivp(lambda t,u: model(t,u)[0],t[[0,-1]],f0,t_eval=t)
并通过再次调用相同的模型函数从结果中得到
z = model(t,result.y.T)[1]
如果模型函数不能自动矢量化,则使用迭代器公式
z = [ model(t_,u_)[1] for t_,u_ in zip(t,result.y.T) ]
z = np.array(z).T # to get the same index ordering as in result.y
如果模型函数只返回一个值列表而不是一对元组,请使用适当的数组切片操作。
,嘿,谢谢你的解决方案。我认为这个解决方案也可能有效:
def model(t,y):
global inner_func
def inner_func(t,y):
#differential Eq System
return dx_dt,dy_dt,z
dx_dt,z = inner_func(t,y)
return dx_dt,dy_dt
t = np.linspace(0,10,100)
f0 = [2,1,0]
result = solve_ivp(model,[np.min(t),np.max(t)],t_eval=t)
z = np.array([inner_func(t[i],result[i,:]) for i in range(len(t))])[:,2]