Processing math: 100%

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
2π01cos(nx)dx
arvoksi saadaan
22n21cos(2nπ)cot(nπ)n.
Jos n on positiivinen kokonaisluku, osoittajaan syntyy 0. Raja-arvotarkastelulla saadaan joko 0 tai 42/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 1cos(nx)=2sin2(nx/2) avulla, ja tulokseksi saadaan 42 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=10000).  Numeerinen integrointi antaa tulokseen 42 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
2π012+cos(x)dx
arvoksi saadaan 2π/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])
23arctan(tan(x/2)3).
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=π 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 12+cos(x) integraalifunktioksi saadaan 23arctan(tan(x/2)3)mod(xπ,2π)x3 missä jatkuvuus on hoidettu aivan oikein, mutta hintana tietenkin on mod-funktion ilmestyminen lausekkeeseen. Integraalin 2π01cos(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:

Anonyymi 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

3i3log(tan(x2)3i)+3i3log(tan(x2)+3i)

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.

Anonyymi kirjoitti...

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