scipy 的 solve_ivp 返回刚性微分方程的错误结果

计算科学 Python scipy 刚性
2021-12-06 00:14:35

我正在使用 scipy 的solve_ivp来求解一个刚性微分方程。我正在使用方法BDF来解决这个问题。我已经使用了 MATLAB 的ode23s,并且在 MATLAB 中得到了正确的结果。但是,solve_ivp 返回错误的结果。问题可能在于计算X我的代码是:

X0=np.array([0,0,0,0,0,0])
swd=0
start =timer()
for k in range(0,100):
    

    tspan=[k*(1e-4), (k+.7)*(1e-4)]
    tspan1=[(k+.7)*(1e-4), (k+1)*(1e-4)]
    
    swm=1
    sol=solve_ivp(conv,tspan,X0,method='BDF')
    TT=sol.t
    YY=sol.y
    
    swm=0
    sol1=solve_ivp(conv,tspan1,YY[:,-1],method='BDF')
    TT1=sol1.t
    YY1=sol1.y
    X0=YY1[:,-1]
    plt.plot(TT,YY[0,:],'r')
    plt.plot(TT,YY[1,:],'b')

end=timer()
print(end-start)

对于定义为的函数

def conv(t,X):
    
    global swd
        if (swd==0) or (swd==1 and X[0]>0):
            
            Xdot = A11.dot(X).reshape(6,1) + B11.dot(U)
            swd=1
        else:
            
            Xdot = A10.dot(X).reshape(6,1) + B10.dot(U)
            swd=0
    return np.squeeze(np.asarray(Xdot))

矩阵A6 6* ,B分别是6 2* 。X6 1*。X之一计算不正确,然后整个图返回不正确。仅正确返回了第一个周期的图。 详细信息:给定的代码是升压转换器。swmswd分别是 MOSFET 和二极管的开关状态。据我调试,二极管在第一个周期正确切换。但是对于第二个循环和其他循环,二极管不会关闭,而是保持打开状态。

0个回答
没有发现任何回复~