我将在下面提供一个可执行的示例,但首先让我带您解决这个问题
我正在使用来自scipy.integrate
的solve_ivp
来解决一个初始值问题(see documentation)。事实上,我必须调用解算器两次,一次向前积分,一次向后积分。(我不得不不必要地深入我的具体问题来解释为什么这是必要的,但请相信我——确实如此!)
sol0 = solve_ivp(rhs,[0,-1e8],y0,rtol=10e-12,atol=10e-12,dense_output=True)
sol1 = solve_ivp(rhs,[0, 1e8],y0,rtol=10e-12,atol=10e-12,dense_output=True)
这里rhs
是初始值问题y(t) = rhs(t,y)
的右侧函数。在我的例子中,y
有六个组件y[0]
到y[5]
y0=y(0)
是初始条件[0,±1e8]
是各自的积分范围,一个在时间上向前,另一个在时间上向后rtol
和atol
是公差
重要的是,您可以看到我标记了dense_output=True
,这意味着解算器不仅返回数值网格上的解,还返回插值函数sol0.sol(t)
和sol1.sol(t)
我现在的主要目标是定义一个分段函数,比如sol(t)
,它取t<0
的值sol0.sol(t)
和t>=0
的值sol1.sol(t)
。因此,主要的问题是:我如何做到这一点?
我认为numpy.piecewise
应该是为我做这件事的首选工具。但我在使用它时遇到了困难,正如你们将在下面看到的,我向你们展示了我迄今为止所做的尝试
下面框中的代码解决了我的示例的初始值问题。大部分代码是rhs
函数的定义,其细节对问题并不重要
import numpy as np
from scipy.integrate import solve_ivp
# aux definitions and constants
sin=np.sin; cos=np.cos; tan=np.tan; sqrt=np.sqrt; pi=np.pi;
c = 299792458
Gm = 5.655090674872875e26
# define right hand side function of initial value problem, y'(t) = rhs(t,y)
def rhs(t,y):
p,e,i,Om,om,f = y
sinf=np.sin(f); cosf=np.cos(f); Q=sqrt(p/Gm); opecf=1+e*cosf;
R = Gm**2/(c**2*p**3)*opecf**2*(3*(e**2 + 1) + 2*e*cosf - 4*e**2*cosf**2)
S = Gm**2/(c**2*p**3)*4*opecf**3*e*sinf
rhs = np.zeros(6)
rhs[0] = 2*sqrt(p**3/Gm)/opecf*S
rhs[1] = Q*(sinf*R + (2*cosf + e*(1 + cosf**2))/opecf*S)
rhs[2] = 0
rhs[3] = 0
rhs[4] = Q/e*(-cosf*R + (2 + e*cosf)/opecf*sinf*S)
rhs[5] = sqrt(Gm/p**3)*opecf**2 + Q/e*(cosf*R - (2 + e*cosf)/opecf*sinf*S)
return rhs
# define initial values, y0
y0=[3.3578528933149297e13,0.8846,2.34921,3.98284,1.15715,0]
# integrate twice from t = 0, once backward in time (sol0) and once forward in time (sol1)
sol0 = solve_ivp(rhs,[0,-1e8],y0,rtol=10e-12,atol=10e-12,dense_output=True)
sol1 = solve_ivp(rhs,[0, 1e8],y0,rtol=10e-12,atol=10e-12,dense_output=True)
从这里可以分别通过sol0.sol
和sol1.sol
处理解决方案函数。作为示例,让我们绘制第四个组件:
from matplotlib import pyplot as plt
t0 = np.linspace(-1,0,500)*1e8
t1 = np.linspace( 0,1,500)*1e8
plt.plot(t0,sol0.sol(t0)[4])
plt.plot(t1,sol1.sol(t1)[4])
plt.title('plot 1')
plt.show()
3.1直接从sol0.sol
和sol1.sol
构建向量值分段函数
def sol(t): return np.piecewise(t,[t<0,t>=0],[sol0.sol,sol1.sol])
t = np.linspace(-1,1,1000)*1e8
print(sol(t))
这导致…/numpy/lib/function_base.py的第628行中出现以下分段错误:
TypeError: NumPy boolean array indexing assignment requires a 0 or 1-dimensional input, input has 2 dimensions
我不确定,但我确实认为这是因为以下原因:在documentation of piecewise中,它提到了第三个论点:
funclistlist of callables, f(x,*args,**kw), or scalars
[...]. It should take a 1d array as input and give an 1d array or a scalar value as output. [...].
我想问题是,我的解决方案有六个部分。因此,在时间网格上计算,输出将是一个2d数组。有人能证实这确实是个问题吗?因为我认为这确实大大限制了piecewise
的实用性
3.2尝试相同的方法,但只针对一个组件(例如第四个组件)
def sol4(t): return np.piecewise(t,[t<0,t>=0],[sol0.sol(t)[4],sol1.sol(t)[4]])
t = np.linspace(-1,1,1000)*1e8
print(sol4(t))
这导致与上述文件相同的第624行出现此错误:
ValueError: NumPy boolean array indexing assignment cannot assign 1000 input values to the 500 output values where the mask is true
与前面的错误相反,不幸的是,到目前为止,我不知道它为什么不起作用
3.3类似尝试,但首先定义了第四个组件的功能
def sol40(t): return sol0.sol(t)[4]
def sol41(t): return sol1.sol(t)[4]
def sol4(t): return np.piecewise(t,[t<0,t>=0],[sol40,sol41])
t = np.linspace(-1,1,1000)
plt.plot(t,sol4(t))
plt.title('plot 2')
plt.show()
现在这不会导致错误,我可以生成一个绘图,但是这个绘图看起来不应该。它应该看起来像上面的图1。同样在这里,我到目前为止还不知道发生了什么
我感谢你的帮助
您可以查看numpy.piecewise源代码。这个函数没有什么特别之处,所以我建议手动执行所有操作
关于你失败的尝试。是,
piecewise
EXPECT函数返回1d数组。您的第二次尝试失败,因为文档说明funclist
参数应该是函数列表或标量,但是您发送了数组列表。与文档相反,它甚至适用于数组,您只需使用与t < 0
和t >= 0
大小相同的数组,如:相关问题 更多 >
编程相关推荐