Sistem gabungan dari 4 persamaan diferensial - Python

Question

Saya sudah mendapatkan sistem gabungan dari 4 persamaan diferensial pada gambar. Saya punya 4 fungsi (xG; yG; gamma; beta) dan turunannya. Mereka semua adalah fungsi dari variabel independen yang sama t.

Saya mencoba menyelesaikannya dengan salep. Masalahnya adalah, untuk melakukannya, saya rasa saya perlu mengungkapkan sistem sedemikian rupa sehingga setiap turunan kedua tidak bergantung pada turunan kedua lainnya. Ini melibatkan sejumlah matematika yang pasti akan membawa saya ke kesalahan di suatu tempat (saya mencoba!).

Tahukah Anda bagaimana saya bisa:

1. Pecahkan sistem persamaan diferensial ini apa adanya?
2. atau dapatkan python untuk mengisolasi turunan kedua untuk saya?

Saya melampirkan kode tes saya

Terima kasih

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

Anda perlu menulis ulang semua turunan orde dua Anda sebagai orde pertama dan menyelesaikan 8 ODE bersama-sama:

Maka Anda membutuhkan kondisi awal untuk semua turunannya, tapi sepertinya Anda sudah punya. FYI, kode Anda tidak berjalan ( line 71: NameError: name 'Xg2' is not defined), silakan periksa.

Juga, untuk informasi lebih lanjut lihat menyelesaikan ODE orde-2 ke-2 secara numerik .

EDIT # 1: Pada langkah pertama, Anda perlu memisahkan sistem persamaan. Meskipun Anda dapat menyelesaikannya secara manual, saya tidak akan merekomendasikan, jadi mari gunakan sympymodul:

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

sekarang kita memiliki semua turunan urutan ke-2 yang dipisahkan. Misalnya, ekspresi untuk beta2adalah

jadi (dan semua turunan urutan kedua lainnya juga) memiliki bentuk

perhatikan bahwa tidak ada ketergantungan pada xgatau yg.

Mari perkenalkan dua variabel baru, bdan k:

kemudian

menjadi

dan sistem lengkap ODE yang harus diselesaikan adalah

Sekarang semua ODE bergantung pada empat variabel yang bukan merupakan turunan dari apapun. Juga, karena xgdan ygdidegenerasi, hanya ada 6 persamaan, bukan 8. Namun, seseorang dapat menulis ulang kedua persamaan ini dengan cara yang sama gammadan betauntuk mendapatkan sistem lengkap dari 8 persamaan, dan mengintegrasikannya bersama.