耦合 OIDENT 求解器中方程的顺序是否重要？

Question

但是，如果我更改公式中定义的方程式的顺序，我的代码也会运行，我的图表也会发生变化。有人能告诉我这是为什么吗？现在我不知道这个系统的正确图表是什么。

def myFunction(r,t):
g = 9.81
L_L = 20 #draught
L_r = 20 #draught
L_d = 4 #ukc
u_s = 0.08
w_d = 60 #width vessel
rho = 1030
b_R = 1.0
b_L = 1.0
b_D = 1.0
A_s = 340*L_L
M_s = 40000*10^3 

w_L = r[0]
w_r = r[1]
u_1 = r[2]
u_2 = r[3]
u_3 = r[4]
u_4 = r[5]
u_d = r[6]
p_1 = r[7]
p_2 = r[8]
p_3 = r[9]
p_4 = r[10]
deltap_L = r[11]
deltap_R = r[12]

u_1 = (L_L*u_s + w_L*u_2)/w_L
u_2 = (w_d*u_d)/w_L
u_d = (w_r*u_3)/w_d
u_3 = (w_L*u_2)/w_r
u_4 = (L_r*u_s + w_r*u_3)/w_r
dwLdt = u_s
dwrdt = - u_s  
du1dt =  - g*((p_2-p_1)/L_L + deltap_L/L_L) - b_L*u_1 
du2dt =  - g*((p_2-p_1)/L_L - deltap_L/L_L) - b_L*u_2
du3dt =  - g*((p_4-p_3)/L_r + deltap_R/L_r) - b_R*u_3
du4dt =  - g*((p_4-p_3)/L_r - deltap_R/L_r) - b_R*u_4
duddt =  - g*((p_3-p_2)/L_d) - b_D *u_d
dp1dt = - u_1
dp2dt = - u_1
dp3dt = + u_4
dp4dt = + u_4
ddeltap_Ldt = -u_1
ddeltap_Rdt = u_4 
#M_s = (g*rho*A_s*((p_1+p_2)/2 + deltap_L/6 - (p_3+p_4)/2 - deltap_R/6))/u_s


return (dwLdt, dwrdt, du1dt, du2dt, du3dt, du4dt, duddt, dp1dt, dp2dt, dp3dt, dp4dt, ddeltap_Ldt, ddeltap_Rdt)
r0 = [10,10,0,0,0,0,0,0,0,0,0,0,0]
t = np.linspace(0,20,10000)
r = odeint(myFunction, r0, t)

w_L = r[:,0]
w_r = r[:,1]
u1 = r[:,2]
u2 = r[:,3]
u3 = r[:,4]
u4 = r[:,5]
ud = r[:,6]
p_1 = r[:,7]
p_2 = r[:,8]
p_3 = r[:,9]
p_4 = r[:,10]
deltap_L = r[:,11]
deltap_R = r[:,12]

plt.figure()
plt.plot(t, u1)
plt.plot(t, u2)
plt.plot(t, u3)
plt.plot(t, u4)
plt.plot(t, ud)
plt.legend(('$u_1$', '$u_2$', '$u_3$', '$u_4$', '$u_d$'))
plt.xlabel('Time ($s$)')
plt.ylabel('Fluid velocity ($m/s$)')

graph first code

如果代码中的方程式是这样的：

    dwLdt = u_s
    dwrdt = - u_s  
    du1dt =  - g*((p_2-p_1)/L_L + deltap_L/L_L) - b_L*u_1 
    du2dt =  - g*((p_2-p_1)/L_L - deltap_L/L_L) - b_L*u_2
    du3dt =  - g*((p_4-p_3)/L_r + deltap_R/L_r) - b_R*u_3
    du4dt =  - g*((p_4-p_3)/L_r - deltap_R/L_r) - b_R*u_4
    duddt =  - g*((p_3-p_2)/L_d) - b_D *u_d
    u_1 = (L_L*u_s + w_L*u_2)/w_L
    u_2 = (w_d*u_d)/w_L
    u_d = (w_r*u_3)/w_d
    u_3 = (w_L*u_2)/w_r
    u_4 = (L_r*u_s + w_r*u_3)/w_r
    dp1dt = - u_1
    dp2dt = - u_1
    dp3dt = + u_4
    dp4dt = + u_4
    ddeltap_Ldt = -u_1
    ddeltap_Rdt = u_4 
    #M_s = (g*rho*A_s*((p_1+p_2)/2 + deltap_L/6 - (p_3+p_4)/2 - deltap_R/6))/u_s

图形是这样的：

graph second mode

为什么会这样？

现在的代码是：

  def myFunction(r,t):
     g = 9.81
     L_L = 20 #draught
     L_r = 20 #draught
     L_d = 4 #ukc
     u_s = 0.08
     w_d = 60 #width vessel
     rho = 1025
     b_R = 1.0
     b_L = 1.0
     b_D = 1.0
     A_s = 340*L_L
     M_s = 40000*10^6

     w_L = r[0]
     w_r = r[1]
     u_1 = r[2]
     u_2 = r[3]
     u_3 = r[4]
     u_4 = r[5]
     u_d = r[6]
     p_1 = r[7]
     p_2 = r[8]
     p_3 = r[9]
     p_4 = r[10]
     deltap_L = r[11]
     deltap_R = r[12]


     u_1 = (L_L*u_s + w_L*u_2)/w_L
     u_2 = (w_d*u_d)/w_L
     u_3 = (w_L*u_2)/w_r
     u_4 = (L_r*u_s + w_r*u_3)/w_r
     u_d = (w_r*u_3)/w_d

     du1dt =  - g*((p_2-p_1)/L_L + deltap_L/L_L) - b_L*u_1 
     du2dt =  - g*((p_2-p_1)/L_L - deltap_L/L_L) - b_L*u_2
     du3dt =  - g*((p_4-p_3)/L_r - deltap_R/L_r) - b_R*u_3
     du4dt =  - g*((p_4-p_3)/L_r + deltap_R/L_r) - b_R*u_4
     duddt =  - g*((p_3-p_2)/L_d) - b_D *u_d




     dp1dt =  - u_1
     dp2dt =  - u_1
     dp3dt = + u_4
     dp4dt = + u_4
     ddeltap_Ldt =  - u_1
     ddeltap_Rdt =  u_4 

     dwLdt = u_s
     dwrdt = - u_s 

     deltap_R = 6*(((p_1+p_2)/2) + (deltap_L/6) - ((p_3+p_4)/2)) 

     return (dwLdt, dwrdt, du1dt, du2dt, du3dt, du4dt, duddt, dp1dt,     dp2dt, dp3dt, dp4dt, ddeltap_Ldt, ddeltap_Rdt)

r0 = [10,10,0,0,0,0,0,0,0,0,0,0,0]
t = np.linspace(0,125,100000)
r = odeint(myFunction, r0, t)

wL = r[:,0]
wr = r[:,1]
u1 = r[:,2]
u2 = r[:,3]
u3 = r[:,4]
u4 = r[:,5]
ud = r[:,6]
p1 = r[:,7]
p2 = r[:,8]
p3 = r[:,9]
p4 = r[:,10]
deltapL = r[:,11]
deltapR = r[:,12]

Answer 1

在代码块中

     u_1 = (L_L*u_s + w_L*u_2)/w_L
     u_2 = (w_d*u_d)/w_L
     u_3 = (w_L*u_2)/w_r
     u_4 = (L_r*u_s + w_r*u_3)/w_r
     u_d = (w_r*u_3)/w_d

您更改了所有 u 值。当然，如果您使用修改或未修改的u 值计算更多值，则函数会有所不同。

最简单的做法是不要重复使用这些变量名。将左侧的u 更改为v，然后检查您是否真的想使用修改后的v 值或u 值。您进行该更改的原始代码将是

     v_1 = (L_L*u_s + w_L*u_2)/w_L
     v_2 = (w_d*u_d)/w_L
     v_3 = (w_L*v_2)/w_r
     v_4 = (L_r*u_s + w_r*v_3)/w_r
     v_d = (w_r*v_3)/w_d

     du1dt =  - g*((p_2-p_1)/L_L + deltap_L/L_L) - b_L*v_1 
     du2dt =  - g*((p_2-p_1)/L_L - deltap_L/L_L) - b_L*v_2
     du3dt =  - g*((p_4-p_3)/L_r - deltap_R/L_r) - b_R*v_3
     du4dt =  - g*((p_4-p_3)/L_r + deltap_R/L_r) - b_R*v_4
     duddt =  - g*((p_3-p_2)/L_d) - b_D *v_d




     dp1dt =  - v_1
     dp2dt =  - v_1
     dp3dt = + v_4
     dp4dt = + v_4
     ddeltap_Ldt =  - v_1
     ddeltap_Rdt =  v_4 

现在对指定类型的任何重新排列都会导致错误，即在某些右侧使用的某些变量以前未定义。

Answer 2

在您所指的第二个代码 sn-p 中，例如u_1 在分配之前。 u_2 ... u_4 也是如此。因此，无论这些值在分配用于计算 du1dt...duddt 之前是什么值