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.