|
很多仿真問(wèn)題的實(shí)質(zhì)都是求解微分方程,，包括常微分方程（ODE）及偏微分方程（PDE）。使用Python的科學(xué)計(jì)算類(lèi)型庫(kù)（numpy,、scipy等）可以方便地對(duì)微分方程進(jìn)行求解,。 對(duì)于簡(jiǎn)單的一階ODE，這里以下落時(shí)的空氣阻力問(wèn)題為例進(jìn)行介紹,。
該問(wèn)題可以由下式進(jìn)行描述：
要進(jìn)行求解,，我們首先將其變換為下式這種形式
變換后得到：
接著我們將上式在python中表達(dá)出來(lái)（系數(shù)任給）
import numpy as npimport matplotlib.pyplot as pltimport scipy as spfrom scipy.integrate import odeintfrom scipy.integrate import solve_ivpdef dvdt(t, v): return 2*v**2 - 3v0 = 0
求解微分方程時(shí)，scipy中有著兩個(gè)常見(jiàn)的求解器：
我們定義時(shí)間為1s，步長(zhǎng)為0.01s,，并使用上述兩種方法分別求解,。 t = np.linspace(0, 1, 100)sol_m1 = odeint(dvdt, y0=v0, t=t, tfirst=True)sol_m2 = solve_ivp(dvdt, t_span=(0,max(t)), y0=[v0], t_eval=t)
將結(jié)果出圖 v_sol_m1 = sol_m1.T[0]v_sol_m2 = sol_m2.y[0]l1,=plt.plot(t, v_sol_m1,'ob')#plot返回的list對(duì)象（list of Line2D）需要解構(gòu)，因此需要在line1和等號(hào)之間加一個(gè)逗號(hào)l2,=plt.plot(t, v_sol_m2, '--r')plt.ylabel('v(t)', fontsize=22)plt.xlabel('t', fontsize=22)plt.legend((l1,l2),['odeint','solve_ivp'])plt.show()
運(yùn)行結(jié)果如下 
|
