Newtonia yleistämään 2

Edellisessä blogikirjoituksessani yleistin Newtonin menetelmän usean yhtälön muodostaman yhtälöryhmän ratkaisemiseen. Yleistämisen tie taitaa kuitenkin olla loputon. Katsotaanpa onnistuisiko neliömatriisin käänteismatriisin hakeminen Newtonin menetelmällä. Tarvittaessa näkökulmaa hieman vaihtaen.

Joitakin pohjatietoja tarvitaan. Tiettyä kokoa olevien neliömatriisien — esimerkiksi \[ A = \begin{pmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{pmatrix} \] — yhteenlasku ja matriisikertolasku (sellaisena kuin se matriisialgebrassa määritellään, ks. esim. https://en.wikipedia.org/wiki/Matrix_multiplication) noudattavat melkein samoja laskusääntöjä kuin reaaliluvut. Poikkeuksena on kertolaskun vaihdannaisuus: matriiseille tulot $AB$ ja $BA$ eivät yleensä ole samoja.

Matriisi $B$ on matriisin $A$ käänteismatriisi, jos $AB = BA = I$, missä $I$ on yksikkömatriisi. Tämän vasemmasta ylänurkasta alkavalla ja oikeaan alanurkkaan päättyvällä laskevalla lävistäjällä on ykkösiä, muualla nollia, esimerkiksi \[ I = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix}.  \] Tällöin merkitään $B = A^{-1}$. Kaikilla matriiseilla ei ole käänteismatriisia.

Matriisien jakolaskua on aina ajateltava käänteismatriisilla kertomisena, kuten reaalilukujen jakolaskukin voidaan ajatella: $c/d = cd^{-1} = d^{-1}c$. Matriiseilla oikealta puolelta ja vasemmalta puolelta jakaminen antavat kuitenkin yleensä eri suuren tuloksen: $CD^{-1} \neq D^{-1}C$, joten jakolaskuja on kaksi erilaista.

Olkoon seuraavassa $A$ neliömatriisi, jolle etsitään käänteismatriisia. Kaikille samankokoisille neliömatriiseille $X$, joilla on käänteismatriisi, voidaan määritellä funktio $F(X) = X^{-1} - A$. Tämän nollakohta, ts. yhtälön $F(X) = O$ ratkaisu on $X = A^{-1}$, ts. $A$:n käänteismatriisi. (Tässä $O$ tarkoittaa pelkistä nollista muodostuvaa neliömatriisia.)

Jotta Newtonin menetelmää voitaisiin soveltaa, on löydettävä funktion $F$ differentiaali, kuten edellisestä blogiartikkelista ilmenee. Tällä kertaa ei kuitenkaan derivointi auta, vaan on vaihdettava näkökulmaa ja pyrittävä esittämään funktion muutos summana lineaarisesta termistä ja korjaustermistä edellisen artikkelin yhtälön \[ f(x) = f(x_0) + f'(x_0)(x - x_0) + \varepsilon(x,x_0)(x - x_0) \] tapaan. Tällöin saadaan seuraavaa: \begin{align*} F(X) &= X^{-1} - A \\ &= (X_0^{-1} - A) - (X_0^{-1}X X_0^{-1} - X_0^{-1}) + (X_0^{-1}X X_0^{-1} - 2X_0^{-1} + X^{-1}) \\ &= (X_0^{-1} - A) - X_0^{-1}(X - X_0)X_0^{-1} + (X_0^{-1} - X^{-1})(X - X_0)X_0^{-1}.  \end{align*} (Algebran taidot ovat tässä tarpeen. Kyseessähän ei ole lausekkeen sievennys vaan sen monimutkaistaminen haluttuun muotoon pääsemiseksi. Lausekkeeseen on lisätty joitakin termejä ja tasapainon säilyttämiseksi vähennetty nämä pois. Täten monimutkaistettua lauseketta on muokattu matriisialgebran säännöillä. Tuloksen tarkistaminen on ehkä helpointa sieventämällä lopputulosta ja päätymällä lähtökohtaan.)

Kolmesta termistä ensimmäinen on funktion arvo pisteessä $X_0$, toinen on lineaarinen termi, ts. $X - X_0$ kerrottuna vakiomatriisilla $X_0^{-1}$ sekä oikealta että vasemmalta. Kolmannessa kaksi ensimmäistä tekijää lähestyvät nollaa, kun $X \to X_0$, ts. se kelpaa korjaustermiksi.

Tällöin kahta ensimmäistä termiä voidaan käyttää funktion $F$ approksimaationa ja yhtälöä $F(X) = O$ likimäärin vastaa yhtälö \[ (X_0^{-1} - A) - X_0^{-1}(X - X_0)X_0^{-1} = O.  \] Kun tämä kerrotaan oikealta ja vasemmalta matriisilla $X_0$, saadaan yhtälö: \[ X_0 - X_0AX_0 = X - X_0.  \] Tästä voidaan ratkaista $X$, jolloin Newtonin menetelmän iteraatiokaavaksi saadaan \[ X_{n+1} = 2X_n - X_nAX_n.  \]

Kovin käyttökelpoinen menettely käänteismatriisin määräämiseen tämä ei kuitenkaan ole, sillä riittävän lähellä olevan aloitusmatriisin löytäminen ei useinkaan ole helppoa. Jos tämä ei ole riittävän lähellä, ei syntyvä iteraatio suppene. Jos matriisi $A$ kuitenkin on lävistäjävaltainen, ts. päälävistäjän alkiot ovat itseisarvoltaan muita riittävästi suurempia, voidaan aloittaa matriisista, jossa lävistäjällä on matriisin $A$ lävistäjäalkioiden käänteisarvot ja muut alkiot ovat nollia. Esimerkiksi matriisin \[ A = \begin{pmatrix} 2 & 1 & 1 \\ 1 & 3 & 1 \\ 1 & 1 & 4 \end{pmatrix} \] tapauksessa alkuarvoksi voidaan ottaa \[ X_0 = \begin{pmatrix} 0.5 & 0 & 0 \\ 0 & 0.33 & 0 \\ 0 & 0 & 0.25 \end{pmatrix}, \] jolloin perättäisiksi käänteismatriisin approksimaatioiksi saadaan \begin{align*} &\left( \begin{array}{ccc} 0.5 & 0 & 0 \\ 0 & 0.33 & 0 \\ 0 & 0 & 0.25 \\ \end{array} \right) \\[3pt] &\left( \begin{array}{ccc} 0.5 & -0.165 & -0.125 \\ -0.165 & 0.3333 & -0.0825 \\ -0.125 & -0.0825 & 0.25 \\ \end{array} \right) \\[3pt] &\left( \begin{array}{ccc} 0.604575 & -0.186467 & -0.1299 \\ -0.186467 & 0.389417 & -0.072402 \\ -0.1299 & -0.072402 & 0.281456 \\ \end{array} \right) \\[3pt] &\left( \begin{array}{ccc} 0.640413 & -0.180814 & -0.121396 \\ -0.180814 & 0.408004 & -0.0618364 \\ -0.121396 & -0.0618364 & 0.291636 \\ \end{array} \right) \\[3pt] &\left( \begin{array}{ccc} 0.646718 & -0.17673 & -0.117862 \\ -0.17673 & 0.411567 & -0.0589877 \\ -0.117862 & -0.0589877 & 0.293982 \\ \end{array} \right) \\[3pt] &\left( \begin{array}{ccc} 0.647058 & -0.176471 & -0.117648 \\ -0.176471 & 0.411764 & -0.058824 \\ -0.117648 & -0.058824 & 0.294117 \\ \end{array} \right) \\[3pt] &\left( \begin{array}{ccc} 0.647059 & -0.176471 & -0.117647 \\ -0.176471 & 0.411765 & -0.0588235 \\ -0.117647 & -0.0588235 & 0.294118 \\ \end{array} \right) \\[3pt] &\left( \begin{array}{ccc} 0.647059 & -0.176471 & -0.117647 \\ -0.176471 & 0.411765 & -0.0588235 \\ -0.117647 & -0.0588235 & 0.294118 \\ \end{array} \right) \end{align*}

Viimeisen approksimaation ja matriisin $A$ tulot antavat todellakin yksikkömatriisin laskentatarkkuudella.

Newtonia yleistämään

Matematiikka tutkii struktuureja. Jos kahdesta asiasta löytyy samanlainen struktuuri, niin toista koskevat tuloksetkin ehkä voidaan siirtää toiseen, mahdollisesti yleistää suppeammasta asiasta laajempaan. Struktuurin tunnistaminen voi kuitenkin edellyttää näkökulman vaihtoa. Katsellaanpa Newtonin iteraatiota:

Yhtälöä $f(x) = 0$ voidaan yrittää ratkaista numeerisesti arvaamalla juurelle $x^*$ jokin likiarvo $x_0$ ja asettamalla funktion $f$ kuvaajalle pisteeseen $(x_0,f(x_0))$ tangentti. Jos $x_0$ on riittävän lähellä juurta $x^*$, tangentin ja x-akselin leikkauspisteestä saadaan tarkempi likiarvo $x_1$. Tangentin yhtälö on $y - f(x_0) = f'(x_0)(x - x_0)$.  Asettamalla tässä $y = 0$ ja ratkaisemalla $x$ saadaan leikkauspisteeksi x-akselin kanssa \[ x_1 = x_0 - \frac{f(x_0)}{f'(x_0)}.  \]

Newtonin menetelmän johto tangentin avulla

Toistamalla menettely arvosta $x_1$ lähtien saadaan uusi likiarvo $x_2$. Kun näin jatketaan, saadaan jono $x_0,\ x_1,\ x_2,\ \dots$, joka varsin yleisillä edellytyksillä suppenee kohden juurta $x^*$. Jono muodostetaan siten iteraatiokaavalla \[ x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}, \quad n = 0,\ 1,\ 2,\ \dots.  \] Menettelyä kutsutaan Newtonin iteraatioksi ja se lienee ensimmäinen numeerinen menetelmä, jonka matematiikan opiskelija oppii.

Ei ole kuitenkaan aivan helppoa nähdä, miten menettely pitäisi yleistää, jos kyseessä onkin yhtälöryhmä. Kahden yhtälön tapauksessa tangenttitasojen kuvittelu voi jotenkin onnistua, mutta isommat yhtälöryhmät ovat jo aika ylivoimaisia mielikuvitukselle.  Näkökulman vaihto kuitenkin auttaa.

 Derivaatta on erotusosamäärän raja-arvo, jolloin lauseke \[ \varepsilon(x,x_0) = \frac{f(x) - f(x_0)}{x - x_0} - f'(x_0) \] lähestyy nollaa, kun $x$ lähestyy arvoa $x_0$. Tämä voidaan saattaa muotoon \[ f(x) = f(x_0) + f'(x_0)(x - x_0) + \varepsilon(x,x_0)(x - x_0).  \] Termiä $f'(x_0)(x - x_0)$ kutsutaan differentiaaliksi, termi $\varepsilon(x,x_0)(x - x_0)$ saa usein nimekseen korjaustermi. Differentiaali on muuttujan $x$ suhteen ensimmäistä astetta, se on ns. lineaarinen termi, joka rajaprosessissa $x \to x_0$ lähestyy nollaa tekijän $x - x_0$ takia. Korjaustermi lähestyy vahvemmin nollaa, sillä sen kumpikin tekijä lähestyy rajaprosessissa nollaa.

Lauseketta $f(x_0) + f'(x_0)(x - x_0)$ voidaan pitää sitä parempana funktion $f(x)$ approksimaationa mitä lähempänä kohtaa $x_0$ ollaan. Voidaankin ajatella, että etsitään tämän nollakohta funktion $f(x)$ nollakohdan sijasta. Aivan sama tämä ei tietenkään ole, mutta ehkä kuitenkin lähempänä. Tällä tavoin päädytään Newtonin menetelmän iteraatiokaavaan. Suppenemisesta ei tämäkään näkökulma vielä kerro, vaan se tarvitsee oman analyysinsa.

Newtonin menetelmän yleistämiseen näkökulma kelpaa. Oleellista on löytää lineaarinen differentiaalitermi ja käyttää sitä funktion approksimoimisessa. Kahden yhtälön tapauksessa ratkaistava yhtälöryhmä on \[ \left\{ \begin{aligned} &f(x,y) = 0, \\ &g(x,y) = 0.  \end{aligned} \right.  \] Kahden muuttujan differentiaalilaskennan mukaan differentiaalitermi on summa funktion osittaisderivaatoista kerrottuna vastaavalla muuttujaerotuksella, funktiolle $f$ siis \[ f_x(x_0,y_0)(x - x_0) + f_y(x_0,y_0)(y - y_0), \] missä $f_x$ ja $f_y$ ovat osittaisderivaatat (ts. funktion derivaatta $x$:n suhteen pitäen $y$:tä vakiona ja derivaatta $y$:n suhteen pitäen $x$:ää vakiona). Vastaavasti funktiolle $g$.

Yhtälöryhmää voidaan nyt approksimoida vastaavasti kuin yhden muuttujan tapauksessa: \[ \left\{ \begin{aligned} &f(x_0,y_0) + f_x(x_0,y_0)(x - x_0) + f_y(x_0,y_0)(y - y_0) = 0, \\ &g(x_0,y_0) + g_x(x_0,y_0)(x - x_0) + g_y(x_0,y_0)(y - y_0) = 0.  \end{aligned} \right.  \] Tämä on lineaarinen yhtälöryhmä tuntemattomina $x$ ja $y$ ja ratkaisu antaa parannetut approksimaatiot $x_1$ ja $y_1$. Toistamalla askelta saadaan generoiduksi jono $(x_0,y_0),\ (x_1,y_1),\ (x_2,y_2),\ \dots$, joka eräin edellytyksin lähestyy yhtälöryhmän ratkaisua $(x^*,y^*)$.

Matriisialgebraan perehtynyt henkilö kirjoittaa lineaarisen yhtälöryhmän muotoon \[ F(X_0) + D_F(X_0)(X - X_0) = O, \] missä \begin{align*} &F(X_0) = \begin{pmatrix} f(x_0,y_0) \\ g(x_0,y_0) \end{pmatrix}, \quad D_F(X_0) = \begin{pmatrix} f_x(x_0,y_0) & f_y(x_0,y_0) \\ g_x(x_0,y_0) & g_y(x_0,y_0) \end{pmatrix}, \\ &X_0 = \begin{pmatrix} x_0 \\ y_0 \end{pmatrix}, \quad X = \begin{pmatrix} x \\ y \end{pmatrix}, \quad O = \begin{pmatrix} 0 \\ 0 \end{pmatrix}.  \end{align*} Ratkaisemalla $X$ saadaan \[ X = X_0 - D_F^{-1}(X_0) F(X_0).  \] Tässä $D_F$ on osittaisderivaattojen muodostama matriisi, joka on Jacobi'n matriisiksi kutsuttu yhden muuttujan funktion derivaatan yleistys. $D_F^{-1}$ on sen käänteismatriisi. Täten on päädytty iteraatiokaavaan \[ X_{n+1} = X_n - D_F^{-1}(X_n) F(X_n), \] jonka analogia yhden muuttujan tapaukseen on ilmeinen. Yleistys useamman yhtälön ryhmiin on suoraviivainen: sama iteraatiokaava pätee, matriisien koko vain on suurempi.

Kuution pyörittelyä

Linkki liikkuvaan kuvaan

xy-tason pistettä $(x,y)$ voi kiertää origon ympäri kulman $w$ verran kertomalla pystyvektoriksi kirjoitetut koordinaatit kiertomatriisilla matriisitulomielessä:
\[
\begin{pmatrix}\cos(w) & -\sin(w) \\ \sin(w) & \phantom{-}\cos(w)\end{pmatrix}
\begin{pmatrix}x \\ y\end{pmatrix}\ .
\]
Kuvion kiertäminen tapahtuu kiertämällä periaatteessa sen jokainen piste, mutta esimerkiksi neliötä kierrettäessä riittää kiertää sen kärkipisteet ja piirtää kierretty neliö näiden avulla.

Mikä olisi kaksiulotteisen kierron kolmiulotteinen vastine? Kiertokeskuksena oleva origo on tällöin korvattava origon kautta kulkevalla kiertoakselilla, jonka ympäri kierretään kulman $w$ verran. Tällöinkin kierto voidaan kuvata matriisilla, joka nyt on kokoa $3 \times 3$.

Matriisi saadaan suhteellisen helposti, jos käytettävissä on riittävästi matriisilaskentaa tukeva symbolinen ohjelma. Jos kiertoakselin yksikön pituinen suuntavektori on $\vec{n} = n_1\vec{i} + n_2\vec{j} + n_3\vec{k}$, muodostetaan matriisi
\[
N = \begin{pmatrix} 0 & -n_3 & n_2 \\ n_3 & 0 & -n_1 \\ -n_2 & n_1 & 0 \end{pmatrix}.
\]
Eksponenttifunktio voidaan määritellä myös matriiseille ja kiertomatriisi saadaan tämän avulla:
\[ Q = \exp(wN),
\] missä $w$ on kiertokulma.

Matriisieksponenttifunktio voidaan määritellä eksponenttifunktion tavanomaisen sarjakehitelmän avulla:
\[
\exp(A) = \sum_{k=0}^\infty \frac{A^k}{k!}.
\]
Tässä potenssit $A^k$ tarkoittavat matriisituloja ja sarjan suppeneminen täytyy luonnollisesti osoittaa. Symboliset ohjelmat eivät kuitenkaan laske matriisieksponenttifunktion arvoja tästä määritelmästä, vaan matriisien ominaisarvoteoria antaa siihen paremmat mahdollisuudet.

Alla oleva lasku näyttää kiertomatriisin laskemisen Mathematicalla tarkkoja arvoja käyttäen. Graafisiin demonstraatiotarkoituksiin on kuitenkin luonnollisempaa käyttää numeerista laskentaa, ts. riittävän tarkkoja likiarvoja.

Mathematica — kuten monet muutkin kehittyneet ohjelmat — tarjoavat mahdollisuuden laatia animaatioita ja säädettäviä demonstraatioita niin ilmiön kuin sitä kuvaavan koodin tutkimiseen: katso kuvakaappausta artikkelin alussa ja linkin takana olevaa liikkuvaa kuvaa.

Matemaattisen tekstin lukijan on syytä olla hereillä ja aina hieman epäluuloinen. Mistä tiedetään, että menettely todella antaa kiertomatriisin? Onko selvää, että kolmiulotteisella kierrolla on aina akseli, kuten edellä ohimennen oletettiin? Jos kuutiota jotenkin pyöräyttelee, niin pääseekö alkuasennosta loppuasentoon aina yhdellä sopivan suuruisella kierrolla sopivan akselin ympäri? Työkaluiksi tarvitaan melkoinen annos matriisilaskentaa. Alaan perehtynyt lukija voi kokeilla taitojaan.