ระบบคู่ของ 4 สมการเชิงอนุพันธ์ - Python

Question

ฉันมีระบบคู่ของสมการเชิงอนุพันธ์ 4 สมการในภาพ ฉันมี 4 ฟังก์ชัน (xG; yG; gamma; beta) และอนุพันธ์ ล้วนเป็นฟังก์ชันของตัวแปรอิสระเดียวกัน t

ฉันกำลังพยายามแก้ปัญหาด้วย odeint ปัญหาคือในการที่จะทำเช่นนั้นฉันคิดว่าฉันจำเป็นต้องแสดงระบบในลักษณะที่อนุพันธ์ของแต่ละวินาทีไม่ได้ขึ้นอยู่กับอนุพันธ์อันดับสองอื่น ๆ สิ่งนี้เกี่ยวข้องกับคณิตศาสตร์จำนวนหนึ่งที่จะพาฉันไปพบข้อผิดพลาดที่ไหนสักแห่ง (ฉันพยายามแล้ว!)

คุณรู้ไหมว่าฉันทำได้อย่างไร:

1. แก้ระบบสมการเชิงอนุพันธ์ตามที่เป็นอยู่?
2. หรือให้ python แยกอนุพันธ์ที่สองให้ฉัน

ฉันกำลังแนบรหัสทดสอบของฉัน

ขอบคุณ

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

คุณต้องเขียนอนุพันธ์ลำดับที่สองของคุณใหม่ทั้งหมดเป็นอนุพันธ์ลำดับแรกและแก้ 8 ODE ด้วยกัน:

จากนั้นคุณต้องมีเงื่อนไขเริ่มต้นสำหรับอนุพันธ์ทั้งหมด แต่ดูเหมือนว่าคุณมีอยู่แล้ว FYI รหัสของคุณไม่ทำงาน ( line 71: NameError: name 'Xg2' is not defined) โปรดตรวจสอบ

นอกจากนี้สำหรับข้อมูลเพิ่มเติมโปรดดูที่การแก้ 2 คำสั่ง ODE ตัวเลข

แก้ไข # 1:ในขั้นตอนแรกคุณต้องแยกระบบสมการออก ในขณะที่คุณสามารถแก้ปัญหาได้ด้วยตนเองฉันไม่แนะนำให้ใช้sympyโมดูล:

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])

ตอนนี้เราได้แยกอนุพันธ์ลำดับที่ 2 ทั้งหมดแล้ว ตัวอย่างเช่นนิพจน์สำหรับbeta2is

ดังนั้น (และอนุพันธ์ลำดับที่ 2 ทั้งหมดด้วย) จึงมีแบบฟอร์ม

ทราบว่าไม่มีการพึ่งพาหรือxgyg

ขอแนะนำตัวแปรใหม่สองตัวbและk:

จากนั้น

จะกลายเป็น

และ ODE ทั้งระบบที่ต้องแก้คือ

ตอนนี้ ODE ทั้งหมดขึ้นอยู่กับตัวแปรสี่ตัวซึ่งไม่ใช่อนุพันธ์ของอะไรเลย นอกจากนี้ตั้งแต่xgและygกำลังถดถอยนอกจากนี้ยังมีเพียง 6 สมการแทน 8. แต่หนึ่งสามารถเขียนทั้งสองสมการในลักษณะเดียวกับgammaและbetaเพื่อให้ได้ระบบเต็มรูปแบบของ 8 สมการและบูรณาการร่วมกัน