Python-oppimateriaali (CHEM-A2600)
scipy.integrate.odeint
Kemian tekniikassa haluamme usein ratkaista (integroida) differentiaaliyhtälöitä numeerisesti. SciPyn integrate-alamoduuli sisältää useita funktioita tätä varten. Tällä kurssilla tutustutaan scipy.integrate.odeint-funktioon.
Otetaan esimerkkinä klassinen esimerkki, eli bakteeripopulaation eksponentiaalinen kasvu. Esimerkkisysteemimme osalta tiedetään, että bakteeripopulaation kasvunopeus dy/dt on suoraan verrannollinen populaation kokoon y ajan hetkellä t:
dy/dt = k * y
Tässä tapauksessa differentiaaliyhtälö on itse asiassa hyvin helppo ratkaista analyyttisesti suoralla integroinnilla (ks. Wikipedia-sivu). Mutta havainnollistetaan tämän suoraviivaisen esimerkin avulla, kuinka differentiaaliyhtälön numeerinen integrointi onnistuu SciPyllä.
Aja alla oleva esimerkki Spyderissä. Huomaa, miten dy/dt on määritelty funktiona f_dy_dt ja miten tämä funktio annetaan odeint-funktion parametriksi. Lisäksi odeint saa parametrinä muuttujan y alkuarvon y_0 ja tutkittavat ajanhetket taulukossa t. Käytännössä siis odeint siis kutsuu funktiota f_dy_dt kaikilla ajan hetkillä t ja antaa parametrinä myös senhetkisen populaation y.
import numpy as np import matplotlib.pyplot as plt import scipy.integrate k = 0.92 # kasvunopeus y_0 = 500 # bakteeripopulaatio alussa def f_dy_dt(y, tx): # Differentiaaliyhtälön määritelmä eksponentiaaliselle kasvulle # Parametri y on populaatio ajan hetkellä t # Huomaa, että aikaparametria tx ei tarvita tässä tapauksessa, mutta # funktion määritelmän täytyy sisältää se, jotta odeint hyväksyy funktion. dy_dt = k * y return dy_dt # Simuloidaan ajanhetket t = [0, 4] h max_t = 4 # tuntia t = np.linspace(0, max_t, max_t * 4 + 1) # pisteet 0.25 h välein # Differentiaaliyhtälön dy_dt numeerinen integrointi # y_0 -> populaation alkuarvo # t -> ajanhetket, joissa populaation määrä lasketaan funktion f_dy_dt avulla # paluuarvo "pop" on taulukko. np.shape(pop) on (len(t), len(y_0)), eli (17,1) pop = scipy.integrate.odeint(f_dy_dt, y_0, t) # Piirretaan kuvaaja teksti = "Bakteeripopulaatio x (k = {}, y_0 = {})".format(k, y_0) plt.plot(t, pop[:, 0], color = 'red', label = teksti) plt.xlim(0, max_t) plt.ylim(0.0, 20000) plt.xlabel('t (h)') plt.ylabel('Bakteeripopulaatio x') plt.legend() plt.show()
Lopputuloksena on seuraavan näköinen kuvaaja:
Huomaa, että odeint-funktion parametri y_0 voi olla myös vektori. Tällöin myös derivaatat laskeva funktio saa parametrikseen vektorin y ja odeint palauttaa taulukon, jossa on yhtä monta saraketta kuin vektorissa y_0 on alkioita.