Toisen kertaluvun differentiaaliyhtälön ratkaisu elementtimenetelmällä

Tarkastellaan toisen kertaluvun differentiaaliyhtälön ratkaisua Galerkinin menetelmällä joka paremmin tunnetaan elementtimenetelmänä. Tarkasteltavana on reuna-arvo-ongelma, jossa yhtälön reunaehtoja esiintyy yhtälön molemmissa päissä ja derivaatat ovat itse asiassa paikkaderivaattoja. Toinen tyyppi olisi alkuarvo-ongelma, jossa derivaatat ovat aikaderivaattoja, mutta siihen tehtävätyyppiin ei tässä nyt mennä.

Valitaan tarkasteltavaksi differentiaaliyhtälöksi \[u''=12x^{2}-36x+18,\qquad u\left(0\right)=u\left(3\right)=0,\qquad(1)\] joka on itse asiassa differentiaaliyhtälön \[ u'''' = 24, \qquad u\left(0\right)=u'\left(0\right)=u\left(3\right)=u'\left(3\right) \] ratkaisu derivoituna kaksi kertaa. Voimme käyttää tätä samaa yhtälö myöhemmissä tutkimuksissa.

Tarkka ratkaisu löytyy integroimalla differentiaaliyhtälöä (eq. 1) kaksi kertaa ja valitsemalla integrointivakiot siten, että reunaehdot toteutuvat. Tarkka ratkaisu siis tiedetään, ja se on \[u\left(x\right)=\left(x-3\right)^{2}x^{2}.\] Tässä esimerkkitapauksessa ratkaisun löytäminen integroimalla oli triviaali toimenpide, mutta yleisesti ottaen osittaisdifferentiaaliyhtälöiden analyyttisiä ratkaisuja ei tunneta ja numeerisia menetelmiä pitää käyttää.

Tarkka ratkaisu on esitetty kuvassa 1. Tarkan ratkaisun voi varmistaa derivoimalla yhtälöä kaksi kertaa sekä tarkastamalla, että reunaehdot toteutuvat.

Figure 1: Differentiaaliyhtälön tarkka ratkaisu.

Haetaan differentiaaliyhtälölle likiratkaisu elementtimenetelmällä. Ratkaisuprosessi lyhyesti kuvattuna menee niin, että kerrotaan differentiaaliyhtälöä testifunktiolla \(v\in\mathcal{V}\) sekä integroidaan alueen yli: \[\int u''\cdot v\,\mathrm{d}x=\int fv\,\mathrm{d}x.\]

On hieman makuasia, pitäisikö ratkaista yhtälöä \(-u''=f\) vai \(u''=f\), monesti mukavuussyistä käytetään miinusmerkkiä. Tämä asia tulee ilmi myöhemmin.

Ratkaisu perustuu osittaisintegrointiin, tai oikeastaan tulon derivointiin. On nimittäin niin, että \[\left(u'v\right)'=u''\cdot v+u'v'.\]

Tästä nyt huomataan, että \[u''\cdot v=\left(u'v\right)'-u'v'.\]

Integroimalla yhtälöä saadaan \[\begin{aligned} \int u''\cdot v\,\mathrm{d}x & =\int\left[\left(u'v\right)'-u'v'\right]\,\mathrm{d}x\\ & =\left[u'v\right]_{x=x_{0}}^{x=x_{1}}-\int u'\cdot v'\,\mathrm{d}x.\end{aligned}\]

Testifunktioavaruus \(\mathcal{V}\) valitaan yleensä siten, että funktiot \(v_{i}\in\mathcal{V}\) häviävät reunoilla, ts. \(v_{i}\left(x_{0}\right)=v_{i}\left(x_{1}\right)=0\), joka johtaa siihen, että sijoitustermi häviää.

Saadaan yhtälön heikko muoto1 eli energiamuoto: määrää \(u\in\mathcal{U}\) siten että kaikille \(v\in\mathcal{V}\) \[\int u'\cdot v'\,\mathrm{d}x=-\int fv\,\mathrm{d}x.\]

Mikäli ratkaistavana olisi ollut yhtälö \(-u''=f\), ei heikossa muodossa olisi tuota miinusmerkkiä. Monesti tämä Poissonin yhtälö kirjallisuudessa esitetäänkin muodossa \[-\Delta u=f,\] missä Laplacen operaattori2 \(\Delta\) määritellään \[\Delta u=\nabla\cdot\nabla u\] ja \(\nabla\) on funktion gradientti.

Lisäksi matemaattisessa kirjallisuudessa heikko muoto esitetään yleensä kanonisessa muodossa: määrää \(u\in\mathcal{U}\) siten, että \[a\left(u,v\right)=l\left(v\right)\qquad\forall v\in\mathcal{V},\] missä \[\begin{aligned} a\left(u,v\right) & =\int u'\cdot v'\,\mathrm{d}x,\\ l\left(v\right) & =\int fv\,\mathrm{d}x.\end{aligned}\]

Elementtimenetelmällä ratkaisu etenee siten, että diskretoidaan ongelma ja valitaan diskretointiin sopiva yriteratkaisu. Esimerkiksi väli \(\Omega=[0,3]\) voidaan diskretoida niin, että kuvataan väli \(\Omega^{\left(1\right)}=[0,2]\) yhdellä funktiolla \(\phi_{1}\left(x\right)\) ja väli \(\Omega^{\left(2\right)}=[2,3]\) toisella.

Tavoitteena on keksiä likiratkaisu \(u_{h}(x)\) joka toteuttaisi yhtälön heikon muodon. Ratkaisussa käytettävät muotofunktiot pitää olla vähintään \(\mathcal{C}^{0}\)-jatkuvia, että heikon muodon vasemman puoleisesta termistä tulee jotakin nollasta poikkeavaa.

Yksinkertaisin valinta on siis paloittain jatkuva lineaarinen funktio. Likiratkaisu alkuperäiseen ongelmaan on \[u_{h}\left(x\right)=\begin{cases} \alpha_{2}\phi_{2}^{\left(1\right)}\left(x\right) & x\in\Omega^{\left(1\right)},\\ \alpha_{2}\phi_{1}^{\left(2\right)}\left(x\right) & x\in\Omega^{\left(2\right)}, \end{cases}\] missä muotofunktiot ovat: \[\begin{alignedat}{1}\phi_{2}^{\left(1\right)} & =\frac{1}{2}x,\\ \phi_{1}^{\left(2\right)} & =3-x. \end{alignedat}\]

Muotofunktiot ovat perinteinen hattufunktio, joka saa arvon 1 solmupisteessä 2, siis kohdassa \(x=2\). Muotofunktiot on esitetty kuvassa fig. 2. Mitä on jäljellä keksiä mikä tuo tuntematon kerroin \(\alpha_{2}\) olisi.

Figure 2: Ratkaisuun käytettävät kantafunktiot.

Kantafunktioiden derivaatat ovat \[\begin{alignedat}{1}\frac{\mathrm{d}}{\mathrm{d}x}\phi_{2}^{\left(1\right)} & =\frac{1}{2},\\ \frac{\mathrm{d}}{\mathrm{d}x}\phi_{1}^{\left(2\right)} & =-1. \end{alignedat}\]

Heikko muoto täytyy integroida kahdessa osassa. Kun approksimaatiot sijoitetaan heikkoon muotoon ja lasketaan integraali, saadaan \[k_{22}^{\left(1\right)}+k_{11}^{\left(2\right)}=-f_{2}^{\left(1\right)}-f_{1}^{\left(2\right)},\] missä \[\begin{aligned} k_{22}^{\left(1\right)} & =\int_{\Omega^{\left(1\right)}}\alpha_{2}\frac{\mathrm{d}}{\mathrm{d}x}\phi_{2}^{\left(1\right)}\cdot\frac{\mathrm{d}}{\mathrm{d}x}\phi_{2}^{\left(1\right)}\,\mathrm{d}x=\frac{1}{2}\alpha_{2,}\\ k_{11}^{\left(2\right)} & =\int_{\Omega^{\left(2\right)}}\alpha_{2}\frac{\mathrm{d}}{\mathrm{d}x}\phi_{1}^{\left(2\right)}\cdot\frac{\mathrm{d}}{\mathrm{d}x}\phi_{1}^{\left(2\right)}\,\mathrm{d}x=\alpha_{2},\\ f_{2}^{\left(1\right)} & =\int_{\Omega^{\left(1\right)}}f\phi_{2}^{\left(1\right)}\,\mathrm{d}x=6,\\ f_{1}^{\left(2\right)} & =\int_{\Omega^{\left(2\right)}}f\phi_{1}^{\left(2\right)}\,\mathrm{d}x=0.\end{aligned}\]

Sijoittelemalla lukuarvoja saadaan \[\begin{aligned} k_{22}^{\left(1\right)} & =\alpha_{2}\int_{0}^{2}\frac{1}{2}\cdot\frac{1}{2}\,\mathrm{d}x=\frac{1}{2}\alpha_{2,}\\ k_{11}^{\left(2\right)} & =\alpha_{2}\int_{0}^{2}\left(-1\right)\cdot\left(-1\right)\,\mathrm{d}x=\alpha_{2},\\ f_{2}^{\left(1\right)} & =\int_{2}^{3}\left(12x^{2}-36x+18\right)\frac{1}{2}x\,\mathrm{d}x=6,\\ f_{1}^{\left(2\right)} & =\int_{2}^{3}\left(12x^{2}-36x+18\right)\left(3-x\right)x\,\mathrm{d}x=0.\end{aligned}\]

Tästä voi nyt ratkaista tuntemattoman kertoimen \(\alpha_{2}\): \[\alpha_{2}=k_{2}^{-1}f_{2}=\frac{2}{3}\cdot6=4.\]

Niinpä \(u\):n likiratkaisu \(u_{h}\) on \[u_{h}\left(x\right)=\begin{cases} 2x & x\in\Omega^{\left(1\right)},\\ 12-4x & x\in\Omega^{\left(2\right)}. \end{cases}\]

Tarkan ratkaisun ja likiratkaisun vertailua on tehty kuvassa 2.

Figure 3: Differentiaaliyhtälön tarkan ratkaisun ja likiratkaisun vertailua.

Elementtimenetelmässä, tai tarkemmin sanottuna Galerkinin menetelmässä3, keskeinen ominaisuus on, että virhe on kohtisuorassa muotofunktioihin nähden, eli \[a\left(u-u_{h},v\right)=0\]

Suoraan integroimalla voi todeta että näin todellakin on: \[\begin{gathered} \int_{0}^{2}\frac{\mathrm{d}}{\mathrm{d}x}\left(\left(x-3\right)^{2}x^{2}-2x\right)\cdot\frac{\mathrm{d}}{\mathrm{d}x}\left(\frac{1}{2}x\right)=0,\\ \int_{2}^{3}\frac{\mathrm{d}}{\mathrm{d}x}\left(\left(x-3\right)^{2}x^{2}-\left(12-4x\right)\right)\cdot\frac{\mathrm{d}}{\mathrm{d}x}\left(3-x\right)=0.\end{gathered}\]

Ratkaisu toimii, mutta siinä on muutamia hankaluuksia numeerista implementointia ajatellen. Eräs mielenkiintoinen kysymys on, kuinka näitä kantafunktioita ja niiden derivaattoja voidaan rakentaa tehokkaasti. Sitä käsitellään tulevissa kirjoituksissa.


  1. https://en.wikipedia.org/wiki/Weak_formulation↩︎

  2. https://en.wikipedia.org/wiki/Laplace_operator↩︎

  3. https://en.wikipedia.org/wiki/Galerkin_method↩︎