python解常微分方程_Python-sympy.dsolve求解常微分方程(组)
這里分別介紹怎么利用sympy.dsolve求解常微分方程和常微分方程組。
#首先利用sympy.dsolve求解單個的常微分方程:
#代碼
from sympy import Function, dsolve, Derivative, symbols
from sympy.abc importt
#sympy.abc表示This module exports all latin and greek letters as Symbols。上行代碼意思是導入字母t,因為t字母在dsolve函數(shù)中要被使用,而dsolve函數(shù)中的x這里不需要導入是因為在下一行會聲明x變量
x = Function('x')
#上行代碼等效于x = symbols('x', cls=Function),意思是聲明x為因變量函數(shù)
result=dsolve(Derivative(x(t),?t, 4) - 22*Derivative(x(t),?t, 2) - 24*x(t))
#result=dsolve(Derivative(x, t, 4) - 22*Derivative(x, t, 2) - 24*x),這樣寫是不可以的
#Derivative(x(t),?t, i)可以寫成diff(x(t),?t, i),或者寫成Derivative(x(t),?t, t, t, t)和diff(x(t),?t, t, t, t)都表示函數(shù)x(t)對t的四階導數(shù),注意如果要使用diff函數(shù),在第一行代碼中要加入from sympy import diff。此外,diff()除了作為函數(shù)使用,還可以作為一個表達式的屬性,如x(t)對t的四階導數(shù)可以表示為x(t).diff(t).diff(t).diff(t).diff(t)
print (result)
#result=Eq(x(t), C1*exp(-t*sqrt(11 + sqrt(145))) + C2*exp(t*sqrt(11 + sqrt(145))) + C3*sin(t*sqrt(-11 + sqrt(145))) + C4*cos(t*sqrt(-11 + sqrt(145))))
#其次利用sympy.dsolve求解常微分方程組:
#代碼
from sympy import Function, dsolve, Derivative, symbols, Eq
#from sympy.abc importt
t=symbols('t')
x, y=symbols('x, y', cls=Function)
#注意這里的x和y之間的逗號可以省略,即可以寫成x, y=symbols('x?y', cls=Function)。這行代碼等效于x=Function('x')和y=Function('y')
eq=(Eq(Derivative(x(t),t, 2), 12*(x(t) + y(t))), Eq(Derivative(y(t),t, 2),?12*x(t) +?10*y(t)))
#注意這里使用了Eq()函數(shù),所以在第一行導入了from sympy import Eq。Eq()函數(shù)內(nèi)的逗號相當于是等于號,等號左邊和右邊分別為微分方程左邊和右邊的表達式,如對于本例的第一個微分方程,除了寫成Eq(Derivative(x(t),t, 2),12*(x(t)+y(t)))外,還可以寫成Eq(Derivative(x(t),t, 2)-12*(x(t)+y(t)), 0),注意這后面的0不能省略。
#同樣地,對于第一個例子求解單個微分方程的情況,如果要使用Eq函數(shù),則應該寫成eq=Eq(Derivative(x(t),?t, 4) - 22*Derivative(x(t),?t, 2) - 24*x(t), 0)和result=dsolve(eq),并且要在第一行代碼導入Eq函數(shù)。
#如果在該例中,我們不想使用Eq()函數(shù),則上行代碼應該寫為eq=(Derivative(x(t), t, 2)-12*(x(t) + y(t)), Derivative(y(t), t, 2)- 12*x(t) - 10*y(t))
#簡而言之就是,如果要使用Eq函數(shù),則需要寫成微分方程等號左右兩邊的表達式,當右邊為0時,也需要寫出來;如果不使用Eq函數(shù),則只需要寫出微分方程等號左邊的表達式,此時等號右邊應為0,且不需要寫出來。
result=dsolve(eq)
print (result)
#result=[Eq(x(t), 12*C1*exp(t*CRootOf(l**4 - 22*l**2 - 24, 0)) + 12*C2*exp(t*CRootOf(l**4 - 22*l**2 - 24, 1)) + 12*C3*exp(t*CRootOf(l**4 - 22*l**2 - 24, 2)) + 12*C4*exp(t*CRootOf(l**4 - 22*l**2 - 24, 3))), Eq(y(t), C1*(-12 + CRootOf(l**4 - 22*l**2 - 24, 0)**2)*exp(t*CRootOf(l**4 - 22*l**2 - 24, 0)) + C2*(-12 + CRootOf(l**4 - 22*l**2 - 24, 1)**2)*exp(t*CRootOf(l**4 - 22*l**2 - 24, 1)) + C3*(-12 + CRootOf(l**4 - 22*l**2 - 24, 2)**2)*exp(t*CRootOf(l**4 - 22*l**2 - 24, 2)) + C4*(-12 + CRootOf(l**4 - 22*l**2 - 24, 3)**2)*exp(t*CRootOf(l**4 - 22*l**2 - 24, 3)))]
#CRootOf(表達式, i)應該就是“表達式=0”的第i+1個根(CRootOf(表達式, i)相當于會形成“表達式=0”的根的列表),這一點我還不確定,可以進一步討論。至于為什么CRootOf中的表達式看著像數(shù)字,但卻不直接計算出來而是寫成式子的形式這一點,我的學Python比較好的同學說是為了語言的形式簡便所以利用了表達式的形式,但是我覺得CRootOf中的表達式的I像是復數(shù)符號,這一點我還不確定,可以進一步討論。
我們注意到其實第一個例子的單個微分方程就是第二個例子的微分方程組通過求導和互相代入得到的,因此第一個例子中求出來的x(t)應該和第二例子求出來的x(t)等價,有興趣的可以驗證這一點(應該會使用到歐拉方程exp(It)=cost+Isint,其中I為復數(shù)符號),這一點我還不確定,可以進一步討論。
總結(jié)
以上是生活随笔為你收集整理的python解常微分方程_Python-sympy.dsolve求解常微分方程(组)的全部內(nèi)容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: java实现rsa欧几里得算法求d_RS
- 下一篇: mysql数据迁移neo4j_neo4j