之前在“使用Python求解微分方程（1）一階ODE”介紹了使用Python的Scipy庫求解一階常系數(shù)微分方程。
今天介紹一個名為Gekko的第三方庫，以及如何使用Gekko求解一階常系數(shù)微分方程。
Gekko介紹
Gekko(壁虎)是一個用于機器學習和動態(tài)優(yōu)化的面向對象的Python第三方庫，它可以與線性、二次、非線性和混合整數(shù)編程（LP、QP、NLP、MILP、MINLP）的大規(guī)模求解器相配合?？梢詫崿F(xiàn)的功能包括參數(shù)回歸、數(shù)據(jù)調節(jié)、實時優(yōu)化、動態(tài)模擬和非線性預測控制。
下面介紹一些Gekko的常規(guī)操作，方便理解后文的方程求解。
使用Gekko時一般需要創(chuàng)建一個模型對象：
from gekko import GEKKOm = GEKKO([server], [name]):
定義常量
c = m.Const(value, [name])#value必須提供且必須是標量定義變量
v = m.Var([value], [lb], [ub], [integer], [name])#lb 和 ub 分別為求解器提供變量下限和上限#integer 是一個布爾值，它為混合整數(shù)求解器指定一個整數(shù)變量
定義方程
m.Equation(equation)#例如m.Equation(3*x == (y**2)/z)
方程均使用隱式求解
求解案例
這里以下述微分方程為例進行求解。
我們令k=0.3，y的初始值y0=5進行計算。
求解代碼為：
import numpy as npfrom gekko import GEKKOimport matplotlib.pyplot as pltm = GEKKO() # create GEKKO modelk = 0.3 # constanty = m.Var(5.0) # create GEKKO variablem.Equation(y.dt()==-k*y) # create GEKKO equationm.time = np.linspace(0,20) # time points# solve ODEm.options.IMODE = 4m.solve(disp=False)# plot resultsplt.plot(m.time,y)plt.xlabel('time')plt.ylabel('y(t)')plt.show()
運行結果為：
接下來可以嘗試修改系數(shù)k，看看結果有何不同。這里令k=0.1、0.5、0.7分別計算試試，求解代碼為：
import numpy as npfrom gekko import GEKKOimport matplotlib.pyplot as pltm = GEKKO() # create GEKKO modelk = m.Param() # constanty = m.Var(5.0) # create GEKKO variablem.Equation(y.dt()==-k*y) # create GEKKO equationm.time = np.linspace(0,20) # time points# solve ODEs and plotm.options.IMODE = 4m.options.TIME_SHIFT=0k.value = 0.1m.solve(disp=False)plt.plot(m.time,y,'r-',linewidth=2,label='k=0.1')k.value = 0.5m.solve(disp=False)plt.plot(m.time,y,'b--',linewidth=2,label='k=0.5')k.value = 0.7m.solve(disp=False)plt.plot(m.time,y,'g:',linewidth=2,label='k=0.7')plt.xlabel('time')plt.ylabel('y(t)')plt.legend()plt.show()
計算結果為：
再嘗試求解一個：
這里令u(t)為一階躍函數(shù)，在t=10時u=1。
求解代碼為：
import numpy as npfrom gekko import GEKKOimport matplotlib.pyplot as pltm = GEKKO() # create GEKKO modelm.time = np.linspace(0,40,401) # time points# create GEKKO parameter (step 0 to 1 at t=10)u_step = np.zeros(401)u_step[100:] = 2.0u = m.Param(value=u_step)y = m.Var(1.0) # create GEKKO variablem.Equation(5 * y.dt()==-y+u) # create GEKKO equation# solve ODEm.options.IMODE = 4m.solve(disp=False)# plot resultsplt.plot(m.time,y,'r-',label='Output (y(t))')plt.plot(m.time,u,'b-',label='Input (u(t))')plt.ylabel('values')plt.xlabel('time')plt.legend(loc='best')plt.show()
計算結果為：
