tiistai 12. huhtikuuta 2016

Mitä symbolilaskentaohjelmalta voi odottaa ja mitä ei, osa 2

Kollega Pekka Alestalo esitti kommentoitavaksi Mathematicalla  lasketun esimerkin, jossa tulos on virheellinen: Integraalin
\[
\int_0^{2\pi} \sqrt{1 - \cos(nx)}\,dx
\]
arvoksi saadaan
\[
\frac{2\sqrt{2}}{n} - \frac{2\sqrt{1 - \cos(2n\pi)}\cot(n\pi)}{n}.
\]
Jos $n$ on positiivinen kokonaisluku, osoittajaan syntyy $0 \cdot \infty$. Raja-arvotarkastelulla saadaan joko $0$ tai $4\sqrt{2}/p$ riippuen siitä, kummalta puolelta $n$ lähestyy kokonaislukuarvoa $p$.

Tuntuu omituiselta, sillä kyseessä on ei-negatiivisen funktion integraali, eikä kokonaislukuparametrissa $n$ pitäisi olla mitään kummallista.  Integraali on mahdollista laskea trigonometrisen kaavan $1 - \cos(nx) = 2\sin^2(nx/2)$ avulla, ja tulokseksi saadaan $4\sqrt{2}$ kokonaislukuparametrin $n$ arvosta riippumatta. Tämä saadaan myös Mathematican avulla antamalla $n$:lle riittävän pieni numeerinen arvo. Liian suuri numeerinen arvo saattaa antaa väärän tuloksen (esimerkiksi $n = 10\,000$).  Numeerinen integrointi antaa tulokseen $4\sqrt{2}$ sopivat likiarvot.

Esimerkki osoittaa, että symboliseen laskentaan on suhtauduttava varovaisesti.  Tulosten järkevyys on vähänkin epäilyttävässä tapauksessa varmistettava esimerkiksi numeerisella laskennalla tai graafisilla tarkasteluilla.  Symboliikka on vaikeata eikä siihen pidä sokeasti luottaa.

Miksi laskenta sitten menee väärin?

Kaupallisen ohjelmiston koodit eivät ole julkisia eikä niistä todennäköisesti olisi kovin helppoa löytääkään syytä. Yksinkertaisempi esimerkki sen sijaan vihjaa mahdolliseen ongelmakohtaan.

Integraalin
\[
\int_0^{2\pi} \frac{1}{2 + \cos(x)}\,dx
\]
arvoksi saadaan $2\pi/\sqrt{3}$ Mathematican komennolla
Integrate[1/(2 + Cos[x]),{x,0,2 Pi}].
Tämä sopii myös yhteen numeerisen integroinnin kanssa, jolloin komento on
NIntegrate[1/(2 + Cos[x]),{x,0,2 Pi}].

Integraalifunktio on (komento Integrate[1/(2 + Cos[x]),x])
\[
\frac{2}{\sqrt{3}}\arctan\left(\frac{\tan(x/2)}{\sqrt{3}}\right).
\]
Kummankin rajan sijoittaminen tähän antaa tulokseksi $0$. Tämän mukaan intgeraalin arvo olisi siis $0$, mikä ilmiselvästi on virheellinen tulos.

Virhe näkyy piirtämällä integraalifunktion kuvaaja:



Funktiolla on hyppyepäjatkuvuus kohdassa $x = \pi$ eikä se siis kelpaakaan sellaisenaan integraalifunktioksi, jonka tulee olla jatkuva. On lisättävä integroimisvakio, erilainen epäjatkuvuuskohdan eri puolilla, jotta saadaan jatkuva funktio.

Alun esimerkissä on samanlainen ongelma. Vaikuttaa siltä, että Mathematican integrointialgoritmi toimii, jos epäjatkuvuuskohdat eivät riipu symbolista eikä niitä ei ole liian paljon.

--------------------
Jatkoa:

Kokeilin esimerkkejä myös TI-Nspire-ohjelmistolla (versio 3.9). Funktion \[ \frac{1}{2+\cos(x)} \] integraalifunktioksi saadaan \[ \frac{2}{\sqrt{3}}\arctan\left(\frac{\tan(x/2)}{\sqrt{3}}\right) - \frac{\text{mod}(x-\pi,2\pi) - x}{\sqrt{3}} \] missä jatkuvuus on hoidettu aivan oikein, mutta hintana tietenkin on mod-funktion ilmestyminen lausekkeeseen. Integraalin \[ \int_0^{2\pi}\sqrt{1-\cos(nx)}\,dx \] laskeminen ei onnistunut, integraalifunktiota ei myöskään löytynyt, ts. vastaukseksi tuli integraali muuttumattomana. Numeerisella $n$:n arvolla vastaukseksi tuli likiarvo.

2 kommenttia:

Nimetön kirjoitti...

Vertailun vuoksi otetaan vielä eräs varsin suosittu avoimen / vapaan lähdekoodin symbolinen kirjasto SymPy:

Erittäin hitaaaaan laskun jälkeen SymPy toteaa ettei se osaa laskea sqrt(1 - cos(nx)):n integraalifunktiota tai määrättyä integraalia.

Yksinkertaisempi 1/(2+cos(x)) sen sijaan saa integraalifunktion

$- \frac{\sqrt{3} i}{3} \log{\left (\tan{\left (\frac{x}{2} \right )} - \sqrt{3} i \right )} + \frac{\sqrt{3} i}{3} \log{\left (\tan{\left (\frac{x}{2} \right )} + \sqrt{3} i \right )}$

joka on kompleksilukuineen ehkäpä hieman vaikeatulkintainen. Sympy osaa kuitenkin piirtää siitä graafin, jossa nähdään samanlainen epäjatkuvuuskohta kuin yllä,

ja näin ollen määrätyksi integraaliksi välillä 0, 2pi SymPy esittää (väärän) 0.

Nimetön kirjoitti...

Niin, ja käytin SymPyn versiota 0.7.5, kirjastosta on tullut myös uudempia versioita, jotka ehkä selviytyvät paremmin.