我得到了图中4个微分方程的耦合系统。我有4个函数(xG;yG;gamma;beta)及其导数。它们都是同一个自变量t的函数
我正试图用odeint解决这个问题。问题是,为了做到这一点,我想我需要用一种方式来表达系统,即每个二阶导数不依赖于其他二阶导数。这涉及到大量的数学知识,肯定会让我在某个地方出错(我试过!)
你知道我怎么能:
我正在附加我的测试代码
谢谢
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)
您需要将所有二阶导数重写为一阶导数,并一起求解8个ODE:
然后你们需要所有导数的初始条件,但看起来你们已经有了。 仅供参考,您的代码没有运行(
line 71: NameError: name 'Xg2' is not defined
),请检查它此外,更多信息见solving the 2nd order ODE numerically
编辑#1: 在第一步,您需要解耦方程组。虽然您可以手动解决它,但我不推荐,所以让我们使用
sympy
模块:现在我们将所有二阶导数解耦。例如,
beta2
的表达式是所以它(以及所有其他二阶导数)的形式是
注意,对
xg
或yg
没有依赖性让我们介绍两个新变量,
b
和k
:然后 变成
要解决的ODE的完整系统是
现在所有的常微分方程都依赖于四个变量,它们不是任何事物的导数。此外,由于
xg
和yg
是退化的,因此也只有6个方程,而不是8个。然而,我们可以用与gamma
和beta
相同的方式重写这两个方程,以获得由8个方程组成的完整系统,并将其集成在一起相关问题 更多 >
编程相关推荐