解常微分方程
解常微分方程問題
例1：假設在平面內有一帶電粒子q，品質為m。空間記憶體在勻強磁場B，方向垂直於平面向內即沿z軸負半軸，以及一個沿y軸負半軸的重力場。帶電粒子從磁場內O點釋放。
則可直接列出粒子的運動方程，將這個方程分解成x和y兩個方向，聯立即可求得該方程組的解。
sympy中的dsolve方法
Python常式
1 #導入 2 from sympy import * 3 import numpy as np 4 import matplotlib.pyplot as plt 5 init_printing() 6 7 ###首先聲明符號x,y,q,m,B,g 8 #參量 9 q,m,B,t,g = symbols('q m B t g',real=True,nonzero=True,positive=True) 10 #函數 11 x,y = symbols('x y',real=True,nonzero=True,positive=True,cls=Function) 12 13 ###再將微分方程表示出來 14 eq1 = Eq(x(t).diff(t,t),-q*B*y(t).diff(t)/m) 15 eq2 = Eq(y(t).diff(t,t),-g+q*B*x(t).diff(t)/m) 16 sol = dsolve([eq1,eq2]) 17 print("微分方程:",sol) #現在列印出sol 18 19 ###這個式子非常冗雜，用trigsimp()方法化簡 20 x = trigsimp(sol[0].rhs) 21 y = trigsimp(sol[1].rhs) 22 print("化簡:",x,y) 23 24 ###將裡面的積分常數算出 25 #定義積分變數，避免報錯,觀察上式輸出式子中有幾個C這裡就定義幾個 26 var('C1 C2 C3 C4') 27 #x(0)=0 28 e1 = Eq(x.subs({t:0}),0) 29 #x'(0)=0,要將subs放在diff後面 30 e2 = Eq(x.diff(t).subs({t:0}),0) 31 #y(0)=0 32 e3 = Eq(y.subs({t:0}),0) 33 #y'(0)=0 34 e4 = Eq(y.diff(t).subs({t:0}),0) 35 l = solve([e1,e2,e3,e4],[C1,C2,C3,C4]) 36 print("積分常數：",l) 37 38 ###將積分常量替換到x和y裡面,我們就得到了解的最終形式 39 x = x.subs(l) 40 y = y.subs(l) 41 print("最終形式：",x,y) 42 43 #作圖 44 ts = np.linspace(0,15,1000) 45 consts = {q:1,B:1,g:9.8,m:1} 46 fx = lambdify(t,x.subs(consts),'numpy') 47 fy = lambdify(t,y.subs(consts),'numpy') 48 plt.plot(fx(ts),fy(ts)) 49 plt.grid() 50 plt.show()
解一階常微分方程：dy/dx=y
Python常式
1 from sympy import * 2 f = symbols('f', cls=Function)#定義函數標識符 3 x = symbols('x')#定義變數 4 eq = Eq(diff(f(x),x,1),f(x))#構造等式，即dy/dx=y 5 #diff(函數,自變數,求導次數) 6 print(dsolve(eq, f(x))) 7 pprint(dsolve(eq, f(x)))#以"pretty"形式列印方程的解
解二階常微分方程：y"+py'+qy=0
Python常式
1 from sympy import * 2 f = symbols('f', cls=Function)#定義函數標識符 3 x,p,q = symbols('x,p,q')#批量定義變數 4 eq = Eq(diff(f(x),x,2)+p*diff(f(x),x,1)+q*f(x),0) 5 #構造方程，即y"+py'+qy=0 6 #diff(函數,自變數,求導次數) 7 print(dsolve(eq, f(x))) 8 pprint(dsolve(eq, f(x)))#以"pretty"形式列印方程的解
常微分方程繪圖
Python常式
1 #繪圖 2 import numpy as np 3 import matplotlib.pyplot as plt 4 from scipy.integrate import odeint 5 6 def diff(y, x): 7 return np.array(1/x) 8 # 上面定義的函數在odeint裡面體現的就是dy/dx =1/x #定義常微分方程 9 x = np.linspace(1, 10, 100) # 給出x範圍，(起始，終值，分割段數) 10 y = odeint(diff, 0, x) # 設函數初值為0，即f(1)=0 11 plt.xlabel('x') 12 plt.ylabel('y') 13 plt.title("y=ln(x)") #定義圖的標題 14 plt.grid()#繪製網格 15 plt.plot(x, y) # 將x,y兩個數組進行繪圖 16 plt.show()
單擺運動：y”+gsin(y)/l=0,g為重力加速度，l為擺長
Python常式
1 from scipy.integrate import odeint 2 import matplotlib.pyplot as plt 3 import numpy as np 4 5 g = 9.8 6 l = 1 7 #重力加速度為9.8m/s,擺長為1m，y"+gsin(y)/l=0,g為重力加速度，l為擺長 8 def diff(d_list, t):#我們可以將一個二階常微分方程分解為含有兩個方程的一階常微分方程組 9 omega, theta = d_list 10 return np.array([-g/l*np.sin(theta), omega]) 11 ''' 12 深入剖析diff函數：diff的左邊代表dω/dt和dθ/dt，由於函數返回的是數組類型，我們 13 可以用這種方式構建一個微分方程組：dθ/dt=ω,dω/dt=-gsin(θ)/l 14 ''' 15 t = np.linspace(0, 10, 1000) 16 result = odeint(diff, [0, 30/180*np.pi], t) 17 # odeint中第二個是初始單擺角度30度，無法採用小角近似 18 plt.xlabel('x') 19 plt.ylabel('y') 20 plt.title("dθ/dt=ω,dω/dt=-gsin(θ)/l") 21 plt.plot(t, result) # 輸出ω和θ隨時變化曲線，即方程解 22 plt.show()
洛倫茲曲線：
以下方程組代表曲線在xyz三個方向上的速度，給定一個初始點，可以畫出相應的洛倫茲曲線。
Python常式
1 import numpy as np 2 from scipy.integrate import odeint 3 from mpl_toolkits.mplot3d import Axes3D 4 import matplotlib.pyplot as plt 5 6 def dmove(Point, t, sets): 7 """ 8 point：present location index 9 sets：super parameters 10 """ 11 p, r, b = sets 12 x, y, z = Point 13 return np.array([p * (y - x), x * (r - z), x * y - b * z]) 14 15 t = np.arange(0, 30, 0.001) 16 P1 = odeint(dmove, (0., 1., 0.), t, args=([10., 28., 3.],)) # 17 ## (0.,1.,0.) is the initial point; args is the set for super parameters 18 P2 = odeint(dmove, (0., 1.01, 0.), t, args=([10., 28., 3.],)) 19 ## slightly change the initial point from y = 1.0 to be y = 1.01 20 21 fig = plt.figure() 22 ax = Axes3D(fig) 23 ax.plot(P1[:, 0], P1[:, 1], P1[:, 2]) 24 ax.plot(P2[:, 0], P2[:, 1], P2[:, 2]) 25 plt.show()
