Hệ thống ghép 4 phương trình vi phân - Python

Question

Tôi đã có hệ thống kết hợp của 4 phương trình vi phân trong hình. Tôi có 4 hàm (xG; yG; gamma; beta) và các dẫn xuất của chúng. Chúng đều là hàm của cùng một biến độc lập t.

Tôi đang cố gắng giải quyết nó bằng thuốc mỡ. Vấn đề là, để làm được như vậy, tôi nghĩ cần phải biểu diễn hệ thống sao cho mỗi đạo hàm cấp hai không phụ thuộc vào các đạo hàm cấp hai khác. Điều này liên quan đến một số lượng toán học chắc chắn sẽ đưa tôi đến một lỗi ở đâu đó (tôi đã thử!).

Bạn có biết làm thế nào tôi có thể:

1. Giải hệ phương trình vi phân này như nó là?
2. hoặc lấy python để cô lập các dẫn xuất thứ hai cho tôi?

Tôi đang đính kèm mã kiểm tra của mình

Cảm ơn

import numpy
import math
from numpy import loadtxt
from pylab import figure,  savefig
import matplotlib.pyplot as plt
# Use ODEINT to solve the differential equations defined by the vector field
from scipy.integrate import odeint



def vectorfield(w, t, p):
    """
    Defines the differential equations for the coupled system.

    Arguments:
        w :  vector of the state variables:
                  w = [Xg, Xg1 Yg, Yg1, Gamma, Gamma1, Beta, Beta1]
        t :  time
        p :  vector of the parameters:
                  p = [m, rAG, Ig,lcavo]
    """
#Xg is position ; Xg1 is the first derivative ; Xg2 is the second derivative (the same for the other functions)
        Xg, Xg1,  Yg, Yg1, Gamma, Gamma1, Beta, Beta1 = w
        Xg2=-(Ig*Gamma2*math.cos(Beta))/(rAG*m*(-math.cos(Gamma)*math.sin(Beta)+math.sin(Gamma)*math.cos(Beta)))
        Yg2=-(Ig*Gamma2*math.sin(Beta))/(rAG*m*(-math.cos(Gamma)*math.sin(Beta)+math.sin(Gamma)*math.cos(Beta)))-9.81
        Gamma2=((Beta2*lcavo*math.sin(Beta))+(Beta1**2*lcavo*math.cos(Beta))+(Xg2)-(Gamma1**2*rAG*math.cos(Gamma)))/(rAG*math.sin(Gamma))
        Beta2=((Yg2)+(Gamma2*rAG*math.cos(Gamma))-(Gamma1**2*rAG*math.sin(Gamma))+(Beta1**2*lcavo*math.sin(Beta)))/(lcavo*math.cos(Beta))
        m, rAG, Ig,lcavo, Xg2,  Yg2, Gamma2, Beta2 = p
    
    
    # Create f = (Xg', Xg1' Yg', Yg1', Gamma', Gamma1', Beta', Beta1'):
    f = [Xg1,
         Xg2,
         Yg1, 
         Yg2, 
         Gamma1, 
         Gamma2, 
         Beta1, 
         Beta2]
         
    return f

    


# Parameter values
m=2.722*10**4
rAG=2.622
Ig=3.582*10**5
lcavo=4
# Initial conditions
Xg = 0.0
Xg1 = 0
Yg = 0.0
Yg1 = 0.0
Gamma=-2.52
Gamma1=0
Beta=4.7
Beta1=0

# ODE solver parameters
abserr = 1.0e-8
relerr = 1.0e-6
stoptime = 5.0
numpoints = 250

#create the time values
t = [stoptime * float(i) / (numpoints - 1) for i in range(numpoints)]
Deltat=t[1]
# Pack up the parameters and initial conditions:
p = [m, rAG, Ig,lcavo, Xg2,  Yg2, Gamma2, Beta2]
w0 = [Xg, Xg1,  Yg, Yg1, Gamma, Gamma1, Beta, Beta1]

# Call the ODE solver.
wsol = odeint(vectorfield, w0, t, args=(p,),
              atol=abserr, rtol=relerr)

Answer 1

Bạn cần viết lại tất cả các dẫn xuất bậc hai của mình dưới dạng các phái sinh bậc nhất và giải 8 ODE cùng nhau:

Sau đó, bạn cần các điều kiện ban đầu cho tất cả các dẫn xuất, nhưng có vẻ như bạn đã có. FYI, mã của bạn không chạy ( line 71: NameError: name 'Xg2' is not defined), vui lòng kiểm tra nó.

Ngoài ra, để biết thêm thông tin, hãy xem giải quyết ODE bậc 2 bằng số .

CHỈNH SỬA # 1: Ở bước đầu tiên, bạn cần giải hệ phương trình. Mặc dù bạn có thể giải quyết nó theo cách thủ công, nhưng tôi không khuyên bạn nên sử dụng sympymô-đun:

import sympy as sm
from sympy import symbols

# define symbols. I assume all the variables are real-valued, this helps the solver. If not, I believe the result will be the same, but just calculated slower
Ig, gamma, gamma1, gamma2, r, m, beta, beta1, beta2, xg2, yg2, g, l = symbols('I_g, gamma, gamma1, gamma2, r, m, beta, beta1, beta2, xg2, yg2, g, l', real = True)

# define left hand sides as expressions
# 2nd deriv of gamma
g2 = (beta2 * l * sm.sin(beta) + beta1**2 *l *sm.cos(beta) + xg2 - gamma1**2 *r * sm.cos(gamma))/(r*sm.sin(gamma))
# 2nd deriv of beta
b2 = (yg2 + gamma2 * r * sm.cos(gamma) - gamma1**2 *r * sm.sin(gamma) + beta1**2 *l *sm.sin(beta))/(l*sm.cos(beta))
# 2nd deriv of xg
x2 = -Ig*gamma2*sm.cos(beta)/(r*m*(-sm.sin(beta)*sm.cos(gamma) + sm.sin(gamma)*sm.cos(beta)))
# 2nd deriv of yg
y2 = -Ig*gamma2*sm.sin(beta)/(r*m*(-sm.sin(beta)*sm.cos(gamma) + sm.sin(gamma)*sm.cos(beta))) - g

# now let's solve the system of four equations to decouple second order derivs
# gamma2 - g2 means "gamma2 - g2 = 0" to the solver. The g2 contains gamma2 by definition
# one could define these equations the other way, but I prefer this form
result = sm.solve([gamma2-g2,beta2-b2,xg2-x2,yg2-y2],
                  # this line tells the solver what variables we want to solve to
                  [gamma2,beta2,xg2,yg2] )
# print the result
# note that it is long and ugly, but you can copy-paste it as python code
for res in result:
    print(res, result[res])

bây giờ chúng ta có tất cả các dẫn xuất bậc 2 được tách rời. Ví dụ, biểu thức cho beta2là

vì vậy nó (và tất cả các dẫn xuất bậc 2 khác cũng vậy) có dạng

lưu ý rằng không có sự phụ thuộc vào xghoặc yg.

Hãy giới thiệu hai biến mới bvà k:

sau đó

trở thành

và toàn bộ hệ thống ODE để giải quyết là

Giờ đây, tất cả các ODE đều phụ thuộc vào bốn biến không phải là dẫn xuất của bất cứ thứ gì. Ngoài ra, vì xgvà yglà suy biến nên cũng chỉ có 6 phương trình thay vì 8. Tuy nhiên, người ta có thể viết lại hai phương trình này theo cách tương tự như gammavà betađể có được hệ phương trình đầy đủ gồm 8 phương trình và tích phân nó với nhau.