计算科学 - 在 scipy 的 solve_ivp 中翻译欧拉码 - 吾爱随笔录

在 scipy 的 solve_ivp 中翻译欧拉码

计算科学 Python scipy

2021-12-25 11:07:13

我的代码基于相似变换X=VZ。我Z通过替换状态空间模型AX+BU和RX+SU变换方程来模拟变换方程的模型，其中 X 被替换为VZ和WZresp 然后我再次应用变换以获得 X 回来。我想使用solve_ivp 执行下面的欧拉代码。虽然我能够构建基本方程，但不能像在欧拉中那样在每个时间步访问解数组以X从Z. X并且Z只是 order 的状态变量矩阵4*1。我的欧拉代码是：

  for i in range(0,500000):
           if (w== 0 and  X[1]> vdon) or (w==1 and X[0]> 0):
                      zdot=inv(V)*A*V*Z+inv(V)*B*U
                      Z(i+1)=Z(i)+ h*zdot
                      w=1
                      X(i+1)=V*Z(i+1)
           else:
                     zdot=inv(W)*R*W*Z+inv(W)*S*U
                     Z(i+1)=Z(i)+ h*zdot
                    w=0
                    X(i+1)=W*Z(i+1)

其中V和是和W的特征向量矩阵分别使用AR

e1,V=LA.eig(A)
e2,W=LA.eig(R)

在 scipy 的 solve_ivp 中，我将函数定义为

def conv(t,Z):
   if (w==0 and X[1]>vdon) or (w==1 and X[0]>0):
      zdot=inv(V)*A*V*Z+inv(V)*B*U
      w=1
   else:
       zdot=inv(W)*R*W*Z+inv(W)*S*U
       w=0
return zdot

此处显示的if条件未按预期工作，仅用于代码理解。我将求解器方程定义为： w=0 #initially sol= solve_ivp(conv, tspan1,X0) aa1=sol.t bb1=sol.y 但我无法在运行时定义X=W*Z或X=V*Z在每个求解器时间步长。

1个回答

无论这是否有意义，您都有一个具有两种操作模式或相位的系统，并且每个相位的相变条件不同，如果系统参数正确，则会启用诸如滞后之类的东西。

所以你有差分系统

def f0(t,Z): return A.dot(Z)+B.dot(U(t))
def f1(t,Z): return R.dot(Z)+S.dot(U(t))

和变化条件

def g0(t,Z): return vdon(t) - V[1].dot(Z)
def g1(t,Z): return W[0].dot(Z)

你的欧拉代码的逻辑然后读作

如果在阶段 0：整合 $Z'=f_0(t,Z)$ 只要 $g_0(t,Z)>0$ . 在该段上，存储 $X(t)=VZ(t)$ . 在 $g_0(t,Z)=0$ 切换到阶段 1。
如果在第一阶段：整合 $Z'=f_1(t,Z)$ 只要 $g_1(t,Z)>0$ . 在该段上，存储 $X(t)=WZ(t)$ . 在 $g_1(t,Z)=0$ 切换到阶段 0。

这个食谱有一个明显的失败，如果两者都 $g_1<0$ 和 $g_2<0$ 那么任何阶段都是不允许的，两个阶段一进入就离开。可以将条件从“虽然 $g_k>0$ “到”直到 $g_k$ 从正到负的变化”，可以使用有方向的事件来实现

g0.terminal=True; g0.direction=-1;
g1.terminal=True; g1.direction=-1;

然后切换机制是

t,z=t0,z0
sw=0
while t<t_final:
    if sw==0:
        sol = solve_ivp(f0,[t,t_final],z,events=g0)
        Z = sol.y.T
        X = V.dot(Z)
    else:
        sol = solve_ivp(f1,[t,t_final],z,events=g1)
        Z = sol.y.T
        X = W.dot(Z)
    # end if
    t,z=sol.t[-1], Z[-1]
    # collect sol.t, Z, X segments in a list to be later concatenated, 
    # see https://stackoverflow.com/a/54796425 (a discontinuous ODE entering sliding mode)
    # or directly append to large lists
#end while

其它你可能感兴趣的问题

上一篇在 (C++) boost::random 和 C++ STL 中，卡方分布比 boost::math 更快吗？下一篇Jos Stam 的“稳定流体”算法与 SIMPLE 和 PISO 算法有什么关系？