分段使用numpy的问题

2024-03-28 12:02:15 发布

您现在位置:Python中文网/ 问答频道 /正文

一,。核心问题

我将在下面提供一个可执行的示例,但首先让我带您解决这个问题

我正在使用来自scipy.integratesolve_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]是各自的积分范围,一个在时间上向前,另一个在时间上向后rtolatol是公差

重要的是,您可以看到我标记了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.solsol1.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()

enter image description here


三,。构建分段函数的尝试失败

3.1直接从sol0.solsol1.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()

enter image description here

现在这不会导致错误,我可以生成一个绘图,但是这个绘图看起来不应该。它应该看起来像上面的图1。同样在这里,我到目前为止还不知道发生了什么

我感谢你的帮助


Tags: 函数outputnppltsolverhssolgm
1条回答
网友
1楼 · 发布于 2024-03-28 12:02:15

您可以查看numpy.piecewise源代码。这个函数没有什么特别之处,所以我建议手动执行所有操作

def sol(t):
    ans = np.empty((6, len(t)))
    ans[:, t<0] = sol0.sol(t[t<0])
    ans[:, t>=0] = sol1.sol(t[t>=0])
    return ans

关于你失败的尝试。是,piecewiseEXPECT函数返回1d数组。您的第二次尝试失败,因为文档说明funclist参数应该是函数列表或标量,但是您发送了数组列表。与文档相反,它甚至适用于数组,您只需使用与t < 0t >= 0大小相同的数组,如:

def sol4(t): return np.piecewise(t,[t<0,t>=0],[sol0.sol(t[t<0])[4],sol1.sol(t[t>=0])[4]])

相关问题 更多 >