perjantai 24. huhtikuuta 2020

Differentiaaliyhtälöiden vaikeus

Differentiaaliyhtälöiden alkeiskurssilla keskitytään usein ratkaisemaan yhtälöiden yksinkertaisia perustyyppejä ja näkemään teoria näiden kautta. Opiskelijalle jää helposti käsitys, että lähes kaikki differentiaaliyhtälöt voidaan ratkaista jollakin algebrallisluonteisella tempulla. Ratkaiseminen tarkoittaa tällöin ns. yleisen ratkaisun löytämistä tavallisten alkeisfunktioiden avulla esitettynä. Tavallisiksi alkeisfunktioiksi kutsutaan potensseja, eksponentti-, logaritmi- ja trigonometrisia funktioita sekä näiden käänteisfunktioita ja niistä peruslaskutoimituksilla tai funktioiden yhdistämisellä saatuja funktioita.

Differentiaaliyhtälön ratkaisun olemassaolo on eri asia. Esimerkiksi varsin yksinkertaiselle differentiaaliyhtälölle $y' = x^2 - y^2$ voidaan piirtää suuntakenttä, koska suoraan yhtälöstä voidaan laskea ratkaisukäyrän derivaatta, ts. sen suunta missä tahansa pisteessä $(x,y)$. Suuntakenttäkuva antaa aiheen uskoa, että jokaisen tason pisteen kautta todellakin kulkee jonkin ratkaisufunktion kuvaaja. Voidaanko funktio lausua alkeisfunktioiden avulla, on täysin eri ongelma. Tietyn pisteen kautta kulkevan (eli tietyn alkuehdon täyttävän) ratkaisufunktion olemassaolon varsinainen todistaminen on hankalampi asia enkä tässä siihen paneudu.


Differentiaaliyhtälön $y' = x^2-y^2$ suuntakenttä
(https://homepages.bluffton.edu/~nesterd/apps/slopefields.html)

Differentiaaliyhtälön $y' = x^2 - y^2$ yleisen ratkaisun lausekkeen löytäminen ei peruskurssitiedoilla onnistu. Hyvät laskentaohjelmat suoriutuvat kuitenkin monista differentiaaliyhtälöistä peruskurssitaitoja paremmin. Esimerkiksi laskentaohjelma Mathematica antaa em. yhtälölle yleisen ratkaisun
\[
y(x)=-\frac{i x^2 \left(-C J_{-\frac{5}{4}}\left(\frac{i x^2}{2}\right)+C
   J_{\frac{3}{4}}\left(\frac{i x^2}{2}\right)-2 J_{-\frac{3}{4}}\left(\frac{i
   x^2}{2}\right)\right)-C J_{-\frac{1}{4}}\left(\frac{i x^2}{2}\right)}{2 x \left(C
   J_{-\frac{1}{4}}\left(\frac{i x^2}{2}\right)+J_{\frac{1}{4}}\left(\frac{i
   x^2}{2}\right)\right)}.
\]
Tässä $C$ on yleisessä ratkaisussa esiintyvä vakio, jonka eri arvoilla saadaan eri pisteiden kautta kulkevat ratkaisukäyrät (ei tosin aivan ongelmitta). Ratkaisussa esiintyy imaginaariyksikkö $i$, vaikka funktion arvot ovatkin reaalisia, kun $x$ on reaalinen. Funktio $J$ eri alaindekseillä on Besselin funktio. Tämä on yksi ns. erikoisfunktioista ja $J_n$ on määritelty eräänä Besselin differentiaaliyhtälön $x^2y'' + xy' + (x^2-n^2)y =0$ ratkaisuna. Tämäkin on sellainen yhtälö, jonka ratkaisua ei voida lausua alkeisfunktioiden avulla. Besselin differentiaaliyhtälö esiintyy monissa fysiikan ilmiöiden mallinnuksissa ja sen seurauksena sen ominaisuudet tunnetaan ja arvot voidaan laskea varsin tarkasti. Laskentaohjelmien kannalta se ei ole sen ihmeellisempi funktio kuin $\sin(x)$ tai $\cos(x)$ (jotka myös voidaan määritellä erään differentiaaliyhtälön ratkaisuina, nimittäin $y'' + y = 0$).

Yleinen ratkaisu kovin monimutkaisen lausekkeen avulla esitettynä ei välttämättä ole käyttökelpoinen.  Sellaistakaan ei läheskään aina voida löytää edes erikoisfunktioiden avulla. Usein onkin luontevampaa etsiä ratkaisua numeerisesti, jolloin tavoitteena on löytää yleisen ratkaisun sijaan vain yksi, sopivan lisäehdon täyttävä ratkaisu.

Tälle hetkelle tyypillinen esimerkki on epidemian kehitystä koskeva SIR-malli. Tämä on kolmen differentiaaliyhtälön ryhmä, jossa on kolme tuntematonta funktiota $S(t)$, $I(t)$ ja $R(t)$, muuttujana aika $t$:
\[
\left\{
\begin{aligned}
S'(t) &=-\beta I(t) S(t), \\ I'(t) &= \beta I(t) S(t) - \gamma I(t), \\ R'(t) &= \gamma I(t).
\end{aligned}
\right.
\]
En puutu tässä yhteydessä siihen, miten malli on muodostettu ja mikä on funktioiden $S$, $I$ ja $R$ merkitys. $\beta$ ja $\gamma$ ovat epidemialle ominaisia vakioita. (Tiivis esitys mallista on esimerkiksi Mikko Rahikan blogissa.  Muutoinkin verkosta löytyy aiheeseen liittyvää materiaalia runsain mitoin.)

Odotin, että Mathematica ei löytäisi mallille yleistä ratkaisua, mutta yllätyksekseni tällainen löytyi. Se on kuitenkin siinä määrin monimutkainen, että käytännöllistä merkitystä ei varmastikaan ole. Jonkinlaisella kritiikilläkin siihen lienee syytä suhtautua.

Numeerisen ratkaisun etsiminen on tarkoituksenmukaisempaa. Tällöin tarvitaan lisäehtoja, esimerkiksi ns. alkuehdot, ts. funktioiden arvot aloitushetkellä. Jos näiksi valitaan $S(0) = 10000$, $I(0) = 1$, $R(0) = 0$ ja vakioille käytetään arvoja $\beta = 0.0007$, $\gamma = 1.5$, saadaan epidemialle tyypilliset ratkaisufunktioiden kuvaajat:

SIR-mallin tyypilliset ratkaisukäyrät

Mathematica antaa ratkaisufunktiot interpoloivina funktioina (InterpolatingFunction). Nämä ovat sopivalla algoritmilla laskettuja approksimaatioita, mutta niitä voidaan Mathematicassa käsitellä funktioina samalla tavoin kuin mitä tahansa muuta funktiota, mm. niitä voidaan derivoida. Jos ratkaisut sijoitetaan alkuperäiseen yhtälöryhmään, tämä ei tarkalleen toteudu, mutta saadaan jonkinlainen kuva approksimaation tarkkuudesta. Alla oleva kuva esittää keskimmäisen yhtälön toteutumisessa olevaa virhettä. Ratkaisujen tarkkuus on käytännössä riittävä, mutta tietty varauksellisuus on paikallaan. Vaikka Mathematican interpoloivia funktioita voidaankin esimerkiksi derivoida, ei derivaattojen tarkkuus ole enää yhtä hyvä.

SIR-mallin keskimmäisen yhtälön toteutumisessa oleva virhe

Mitä tästä pitäisi oppia? Laskentaohjelma on hyvä työkalu, mutta ei pidä odottaa, että se ratkaisee kaikki matemaattiset ongelmat vaivatta. Käyttäjän tulee ymmärtää, mitä on tekemässä.