76315A ATKII - NUMEERINEN MALLINTAMINEN

SISÄLLYSLUETTELO:

Mathematica ja vinkkejä sen kanssa elämiseen.....3

Funktiot............................................................................16

Yhtälöt...............................................................................19

Taulukoiden muodostaminen ja interpolointi	........24

Mittaustulosten käsittelyä............................................31

Raja-arvot ja sarjat.........................................................34

Derivointi ja integrointi.................................................36

Differentiaaliyhtälöitä Mathematicalla.......................42

Mekaniikkaa.....................................................................51

Sähköoppia......................................................................52

Populaatioiden kilpailu.................................................55

Ohjelmointi Mathematica'lla.........................................61

Schrödingerin yhtälö.....................................................76

Vektoreita ja matriiseja..................................................86

Mathematica ja vinkkejä sen kanssa elämiseen

Lausekkeet, symbolit ja komennot

Unix-systeemissä: Mathematica käynnistyy komennolla mathematica,  jolloin avautuu erillinen X-ikkuna,  graaffinen käyttäjäliittymä (notebook),  Mathematica' n komentoihin. Istunto lopetetaan valitsemalla FILE valikosta Quit. Mathematica on myös mahdollista käynnistää komennolla math, jolloin erillistä ikkunaa ei avaudu.

Mathematica' n komennot kirjoitetaan ikkunaan ja ne suoritetaan painamalla  SHIFT ja ENTER nappuloita yhtä aikaa.  Mathematica numeroi suorittamansa käskyt  In[n] ja niihin annetut vastaukset  Out[n] juoksevin numeroin siten, että käskyllä ja vastauksella on sama numero.

1 + 1

Aiempiin tuloksiin voi viitata %n-merkinnällä, missä n on tulosterivin numero. Esimerkiksi

%1 + 4

Toinen mahdollisuus on käyttää %-merkkiä viittaamaan viimeisimpään tulokseen, %% - merkkejä toiseksi viimeiseen tulokseen ja niin edelleen.

Lausekkeet ja symbolit

Mathematicassa käytetään kirjaimia ja numeroita symbolien nimien esittämiseen, kuten matematiikassakin. Symbolien nimissä isot ja pienet kirjaimet tarkoittavat eri asiaa ja vain englannin kielen kirjaimet kelpaavat. Nimellä ei ole pituusrajoitusta, mutta se ei saa alkaa numerolla. Siten x2 on hyväksyttävä symboli, mutta 2x ei ole, vaan se tulkitaan kertolaskuksi 2*x.

Yksinkertainen tapa antaa symbolille arvo on sijoituslause.
Sijoituslauseita on kahta tyyppiä
                     symboli=lauseke	(suora asetus)
                     symboli:=lauseke	(viivästetty asetus)
Operaattori  = laskee lausekkeen niin pitkälle kuin mahdollista, ja määrittelee saadun tuloksen symboliksi,  ns. suora asetus.  Operaattori  := määrittelee lausekkeen sellaisenaan symboliksi, ja sen arvo lasketaan aina uudelleen, kun symbolia käytetään (viivästetty asetus).

Esimerkiksi

a = 5

b := a + 1

b

a = 6

b

d = 3

c = a + b + d

c

d = 1

c

Ensimmäisessä osassa  a:n arvon muuttaminen muuttaa b:n arvon. Jälkimmäisessä tapauksessa näin ei tapahdu.  Symbolin arvon voi tulostaa aina kirjoittamalla sen nimi tai käyttämällä Print-komentoa. 
Esimerkiksi

Print[a]

Symbolin arvon voi hävittää merkinnällä
        x=.
tai
        Clear[x]
Sen jälkeen symbolia x voi käyttää esimerkiksi muuttujana.

Mathematicassa käytetään seuraavia merkintöjä peruslaskutoimituksille.
	Kertolasku  a*b,  a b     (HUOM! välilyönti)
	Jakolasku a/b 
	Yhteen- ja vähennyslasku  a+b, a-b
	Potenssi  a^b
Kertolasku voidaan merkitä siis kahdella tavalla. a*b ja   a b. Miksi tarvitaan välilyönti? Yksinkertaisesti siksi, että muuten Mathematica ei tiedä, tarkoitammeko oliota nimeltä ab vai a:n ja  b:n tuloa.

Komennot

Komennot Mathematicassa ovat aina muotoa 
              Komento[syötteet]
Mathematica on symbolisen laskennan ohjelmisto. Symbolinen laskenta tarkoittaa, että yhtälöitä voi esimerkiksi sieventää tai niistä voi yrittää eliminoida muuttujia. Symboliset ratkaisut ovat aina tarkkoja. Tutki seuraavaa esimerkkiä, jossa ratkaistaan toisen asteen yhtälö 2x^2+3x+3=4.

Clear[x]
Solve[2  x^2 + 3  x + 3 == 4, x]

Toisen asteen yhtälöä kirjoitettaessa käytettiin == merkkiä yhtäsuuruuden ilmoittamiseen. Yhtälön ratkaisut esitetään sääntöinä, jotka on sijoitettu ryhmäsulkujen ({}) sisään. Ryhmäsulkujen avulla muodostettuja olioita kutsutaan listoiksi.

Listat

Lista on joukko alkioita, jotka on erotettu toisistaan pilkuilla ja kirjoitettu aaltosulkujen sisään.
Esimerkiksi

lista1  =  {9, 8, 7, 6, 5}

matriisi = {{K, L, M}, {N, O, P}, {Q, R, S}}

Symboli lista1 on lista, jonka alkiot ovat lukuja. Matematiikassa sillä voidaan esittää pystyvektoria. Symboli matriisi viittaa listaan, jonka alkiot ovat listoja. Tässä tapauksessa kyseessä voisi olla 3x3 matriisi, jonka alkiot on lueteltu vaakariveittäin.  
Listojen alkioihin viitataan kaksinkertaisilla hakasuluilla.

lista1[[3]]

matriisi[[2]]

matriisi[[1, 3]]

Merkintä matriisi[[1,3]] viittaa siis ensimmäisen alkion kolmanteen jäseneen.

Säännöt

Solve-komennon tulos saatiin siis sääntöinä

Solve[2  x^2 + 3  x + 3 == 4, x]

Säännöissä korvataan nuolen vasemmalla puolella olevat oliot x nuolen oikealla puolella olevilla olioilla silloin, kun sääntöä halutaan käyttää.  Esimerkiksi halutaan laskea lausekkeen x^3+5 arvo ensimmäistä sääntöä käyttäen.

x^3 + 5 /. %[[1]]

Lausekkeen ja säännön väliin kirjoitetaan sijoitusoperaattori
                 /.
Symboli x on lausekkeessa x^3+5 muuttuja, jolle sijoitetaan arvo edellisen listan ensimmäistä sääntöä käyttäen. 

Lausekkeet

x1  +  x2  +  x3  /.  
      x1 -> 2

A x1  +  B x2  +  
        C x3 /. {B -> A, 
        C -> A}

ovat myös esimerkkejä lausekkeista, joissa esiintyy sääntö tai lista säännöistä. Ensimmäisessä esimerkissä korvataan x1 arvolla 2 ja toisessa esimerkissä {B->A,C->A} on lista, jossa on kaksi sääntöä, ``korvaa B A:lla, korvaa C A:lla''.

Tehtäviä

1. Laske 3sin(4.0). Lisää tulokseen 3. Jaa tulos, joka sinulla oli ennen 3:n summaamista, 2:lla.
2. Kokeile, mitä tekevät  a==a, 2==1, a=2.
3 .Mieti,miten voit sijoittaa säännön x->Sqrt[2]+5 mukaisen arvon muuttujalle x pysyväksi arvoksi.
4. Kokeile suoran ja epäsuoran sijoituksen eroja sijoittamalla x=a, suora=2x, viive:=2x ja x=(k+1). Mitkä ovat nyt muuttujien suora ja viive saamat arvot?
5. Ratkaise yhtälö (y^3)+4y-9=0.

Luonnonvakiot ja kompleksiluvut

Mathematica tuntee joukon  luonnonvakioita  ja erikoisempia matemaattisia olioita, kuten Neperin luvun e ja π:n. Ohessa lista muutamasta yleisimmästä:

E

Pi

Infinity

I

Kompleksiluvut  esitetään muodossa a + b I , missä a ja b ovat reaalilukuja. 
Esimerkiksi lasketaan kahden kompleksiluvun tulo

a = 2 + 3  I
b = 3 - I
a b

Lausekkeen numeerinen arvo

Komennolla
             N[lauseke] tai lauseke//N
lasketaan lausekkeen arvo desimaalilukuna, niiltä osin kuin se on mahdollista. 
Esimerkiksi

Sqrt[2] + 5/(Pi^2  +  10)

N[%]

E^3/5 // N

Lausekkeiden sieventäminen ja muokkaus

Lausekkeiden sieventämiseen ja muokkaukseen on runsaasti komentoja.
Seuraavassa luettelo muutamista tärkeimmistä:
Expand      	suorita kertolaskut ja potenssiin korottamiset
Factor      	etsi yhteiset tekijät
Together   	 kirjoita lauseke yhteistä nimittäjää käyttäen
Apart       	hajoita termit yksinkertaisia nimittäjiä käyttäen
Cancel      	supista yhteiset tekijät osoittajasta ja nimittäjästä
Simplify    	kokeile algebrallisia muunnoksia ja etsi niistä yksinkertaisin 
Numerator   	osoittaja
Denominator 	nimittäjä

Lauseke, johon komento kohdistuu kirjoitetaan hakasulkeisiin tai käytetään merkintää //. Esimerkiksi
               Expand[lauseke]
tai
               lauseke//Expand

e = (x - 1)^2 (2 + 
              x)/((1 + x) (x - 3)^2)

\(\(\(	\)\(f = Expand[e]\)\)\)

\(\(\(	\)\(Together[f]\)\)\)

\(\(\(	\)\(Apart[%]\)\)\)

\(\(\(	\)\(Factor[%]\)\)\)

\(\(\(	\)\(Numerator[%]\)\)\)

\(\(\(	\)\(Expand[%]\)\)\)

\(\(\(	\)\(x/%\)\)\)

\(\(\(	\)\(Denominator[%]\)\)\)

Kuvat

Yksi Mathematican vahvoja ominaisuuksia on sen sisäänrakennettu kyky tuottaa kuvia. Kaksi- ja kolmiulotteisten kuvien tekemiseen löytyy monia komentoja. Näistä yksinkertaisin, mutta hyvin käyttökelpoinen, on Plot. 
Esimerkkinä harjoitus suhteellisuusteorian kurssista , jossa piti tutkia suureen

c =. 
v =.

1/Sqrt[1 - (v/c)^2]

käyttäytymistä suhteellisen nopeuden v/c funktiona. Tämän teemme piirtämällä kuvan.

Plot[1/Sqrt[1 - a^2], {a, 0, 1}]

Muita kuvanpiirtokomentoja ovat Plot3D, joka piirtää kolmiulotteisen kuvan...

Plot3D[
      Sin[x y] Cos[x] Exp[y], {x, -2, 2}, {y, -2, 2}]

...sekä ContourPlot, joka piirtää tasa-arvokäyrän

ContourPlot[
      Sin[x y] Cos[x] Exp[y], {x, -2, 2}, {y, -2, 2}]

Apua

Apua lyötyy HELP valikosta. Koko Mathematica kirja on käytössä hypertekstien avulla.

Optiot

Monet mathematican komennoista hyväksyvät erilaisia optioita, jotka määräävät esimerkiksi numeerisen tarkkuuden, kuvioiden piirtovärin  ja -tarkkuuden sekä muita enemmän tai vähemmän hyödyllisiä asioita. Jonkin komennon optiot saa selville komennolla Options}.

Esimerkiksi Plot3D valitsimet listautuvat alla olevalla komennolla

Options[Plot3D]

Tehtäviä

6. Laske sin(π) ja cos(π).
7. Laske (Sqrt[-1]), 	i^2, 	i^3 ,   	i^4.
8. Ratkaise yhtälö (z^4)= -1.
9. Piirrä käyrä y=x välillä [-10,2].
10. Hae näkyville Sin-komennon syntaksi.
11. Kokeile mitä ?? tekee! Käytä esimerkkinä vaikka Plot-komentoa.
12. Tutki, osaako Mathematica antaa tietoja jonkin komennon optioiden käytöstä ?-komennon avulla.

Ohjelmat

Mathematica-ohjelmia käsitellään tarkemmin luvussa Ohjelmointi. Mainittakoon kuitenkin, että Mathematica-komentoja voi kirjoittaa tiedostoon, jonka Mathematica sitten suorittaa.

Funktiot

Valmiit funktiot

Mathematicassa on kaikki matamatiikan standardifunktiot.
Abs[x], Sqrt[x], Exp[x], Log[x], Log[a,x]

Sin[x], Cos[x], Tan[x], Cot[x], Sec[x], Csc[x]
ArcSin[x], ArcCos[x], ArcTan[x], ArcCot[x], ArcSec[x], ArcCsc[x] 

Sinh[x], Cosh[x], Tanh,[x] Coth[x], ...
ArcSinh[x], ArcCosh[x], ArcTanh[x], ArcCoth[x], ...

n!, n!!, Binomial[n,m]

Max[k,l,m,...], Min[k,l,m,...]
Round[x], Floor[x], Ceiling[x], Sign[x]

Re[z], Im[z], Conjugate[z], Abs[z], Arg[z]

Random[], SeedRandom[]
...
Käytettävissä on myös suuri joukko erikoisfunktioita
LegendreP[n,x]
SphericalHarmonicY[l,m,theta,phi]
ChebyshevT[n,x]
HermiteH[n,x]
LaguerreL[n,x]
JacobiP[n,a,b,x]
BesselI[n,z]
BesselJ[n,z]
Zeta[s,a]

Esimerkkinä piirretään Chebyshev'n polynomi (Subscript[T, 1])

Plot[ChebyshevT[1, x], {x, -1, 1}]

Omat funktiot

Mathematicassa määritellään funktio muodossa
        nimi[x_,y_,z_,...] := lauseke
tai
        nimi[x_,y_,z_,...] = lauseke
Sijoitusperaattorit := ja = toimivat samoin kuin aikaisemmin kuvattiin. Edellinen on viivästetty sijoitus ja jälkimmäinen suora sijoitus. Funktion nimi on  nimi ja argumentit  x, y, z,... luetellaan alaviivalla varustettuna. Alaviivaa käytetään muuttujien nimeämisessä, mutta sitä ei käytetä, kun lauseketta muodostetaan argumenttien avulla.
Esimerkkinä määritellään funktio 
		f(x,y)=(x^3)+4 sin(π/x)+y 
ja lasketaan sen arvo pisteessä (3,1).

Clear[f]

f[x_, y_] := x^3  +  4 Sin[Pi/x]  +  y

f[3, 1]

Funktion arvon jollakin muutujan arvolla saa laskettua sijoittamalla haluttu arvo argumentiksi

Tehtäviä

13. Määrittele funktio f(x)=x^2+2x. Laske funktion arvo kun x=3.
14. Etsi edellä määritellyn funktion nollakohdat.
15 . Määrittele f(x)=(x^4)-3x+(e^(-x))  Mathematica'ssa ja piirrä sen kuvaaja välillä [-3,3]. 
16. Määrittele g(x)=(x^2)-3x ja ratkaise yhtälö g(x)=4.
17. Määrittele funktio z(x,y)=sin(x y) cos(x y) ja piirrä se 3D-kuvana ja tasa-arvokäyrinä.

Kertaustehtäviä

K1. Määrittele funktio f(x)=x^5-3x^4+12x^2+3.
K2. Piirrä ko. funktion kuvaaja välillä [-5,3]
K3. Ratkaise funktion nollakohdat (Älä säikähdä omituisen näköistä vastausta. Se muuttuu mukavamman näköiseksi komennolla N[%]. Kts. helppi-tiedostosta lisätietoja Solvesta, jos haluat tietää mitä vastauksen merkinnät tarkoittivat). Tarkista lisäksi toteuttaako reaalinen ratkaisu yhtälön f(x)=0 (tästä huomaa, että komento N[%] todella ottaa numeerisen arvon nollakohdista).

Yhtälöt

Polynomiyhtälöt

Mathematica osaa ratkaista Solve komennolla polynomiyhtälöitä. Tarkat ratkaisut löytyvät kaikille astetta k≤4 oleville polynomiyhtälöille, mutta joskus ratkaisut ovat hyvin monimutkaisen näköisiä.  Numeerisesti samat yhtälöt voidaan ratkaista komennolla NSolve.   Solve-komennossa annetaan ratkaistava yhtälö  (tai joukko yhtälöitä pilkulla erotettuna aaltosuluissa), ja muuttuja, (tai joukko muuttujia pilkulla erotettuna aaltosuluissa) jonka suhteen ratkaisu halutaan. 
Esimerkiksi ratkaise yhtälöryhmä
		x  +   y = 0
		x - 2 y = 3

Solve[{x + y == 0, x - 2 y == 3}, {x, y}]

NSolve[{x + y == 0, x - 2 y == 3}, {x, y}]

Reduce-komento ratkaisee samat yhtälöt kuin Solve, mutta esittää tuloksen mahdollisimman yleisessä muodossa. 
Esimerkiksi ratkaistaan yleinen toisen asteen yhtälö a(x^2) + bx + c = 0.

a =. 
b =. 
c =.

Solve[a x^2  +  b x  +  c == 0, x]

Reduce[a x^2  +  b x  +  c == 0, x]

Tuloksen esittämisessä Reduce käyttää loogisia operaattoreita, jotka ovat samoja kuin C-kielessä.
         !      ei
         ||     tai
         &&     ja

Algebrallisen yhtälön numeerinen ratkaiseminen

Yhtälöt, jotka eivät ole polynomiyhtälön muotoa, joutuu yleensä ratkaisemaan numeerisesti.  NSolve-komentoa käytetään polynomiyhtälöiden numeeriseen ratkaisemiseen, kuten edellä esitettiin, mutta  FindRoot-komennolla voidaan ratkaista yleisempiä yhtälöitä. 
Muotoa
         FindRoot[vasen==oikea,{x,x0}]
käytetään silloin, kun yhtälön voi derivoida symbolisessa muodossa (Newtonin menetelmä). Symboli x0 on ensimmäinen arvaus ratkaisulle. Sen keksimisessä on suureksi avuksi, kun piirtää Plot-komennolla sekä vasemman että oikean puolen kuvaajan ja etsii niiden leikkauspisteitä. FindRoot löytää yleensä arvausta lähinnä olevan ratkaisun iteroimalla eli parantamalla arvausta kierros kierrokselta. 
Ratkaistaan esimerkkinä yhtälö tan/x)=x.

FindRoot[Tan[x] == x, {x, 0.5}]

Tarkka ratkaisu on x=0, joten tulos tuntuu epätarkalta, mutta oikean ja vasemman puolen erotus on kuitenkin hyvin pieni

Tan[x] - x /. %

Tarkkuutta voidaan lisätä kasvattamalla parametrin AccuracyGoal arvoa ja lisäämällä iteraatioiden lukumäärää MaxIterations parametrin avulla.

FindRoot[Tan[x] == x, {x, 0.5}, AccuracyGoal -> 15, 
      MaxIterations -> 50]

Tan[x] - x /. %

Toinen mahdollisuus on etsiä ratkaisua käyttämällä kahta lähtöpistettä ja niiden kautta piirrettyä suoraa (sekanttimenetelmä). Uusi parempi arvaus saadaan etsimällä suoran nollakohta ja käyttämällä sitä seuraavalla iterointikierroksella toisen pisteenä. Kaikki tämä tapahtuu tietysti FindRoot-komennon sisällä ja ulos tulee vain ratkaisu. Näin täytyy menetellä silloin, kun yhtälön derivaattaa ei ole mahdollista laskea symbolisessa muodossa.
         FindRoot[vasen==oikea, {x, x1, x2}]
Ratkaistaan sama yhtälö kuin edellä tan(x)=x.

FindRoot[Tan[x] == x, {x, -0.1, 0.5}]

Lineaarinen yhtälöryhmä

Lineaariset yhtälöryhmät muodostavat oman erikoisryhmänsä, joiden ratkaisemisessa kannattaa käyttää LinearSolve-komentoa. 
         LinearSolve[m,b]
Ratkaisuna saadaan vektori x, joka on matriisiyhtälön 
		m  x = b
ratkaisu. Lineaarisen yhtälön kertoimet täytyy siis  antaa matriisimuodossa ja yhtälön oikea puoli vektorimuodossa. Esimerkiksi ratkaistaan lineaarinen yhtälöryhmä
        3 x1 + x2 = 2
       x2 + 2 x3 = 2
        x1 + x3   = 1
Yhtälön kerroinmatriisi syötetään Mathematicalle muodossa

m = {{3, 1, 0}, {0, 1, 2}, {1, 0, 1}}
b = {2, 2, 1}

MatrixForm[m]

MatrixForm[b]

LinearSolve[m, b]

Tehtäviä

18. Ratkaise yleinen kolmannen asteen yhtälö ax^3+bx^2+cx+d=0 käyttämällä reduce-komentoa.
19. Etsi yhtälön cos(x)=x ratkaisu (Vinkki: piirrä ensin kuvaaja cos(x)-x vaikkapa välillä [-6,6], niin näet mihin nollakohta asettuu.).
20. Ratkaise seuraava yhtälöryhmä
	3x+y=2 
	y+2z=2 
	x+z=1
symbolisesti ja numeerisesti.
 
21. Etsi yhtälöparille
	x^2+3 x-1= 0
	3 x^2+5 x+1 = 0
symboliset ja numeeriset ratkaisut.

22. Ratkaise yhtälöryhmä
	x^2+5 x-4 y=0
	2 x-y=3

23. Ratkaise lineaarinen yhtälöryhmä
	3 (Subscript[x, 1])+(Subscript[x, 2])+4 (Subscript[x, 3])=2
	5(Subscript[x, 1])+(Subscript[x, 2]) +2(Subscript[x, 3]) =3
	(Subscript[x, 1])-3 (Subscript[x, 2]) -2(Subscript[x, 3]) =1

Taulukoiden muodostaminen ja interpolointi

Fysiikassa ja kemiassa tulokset esitetään usein taulukon muodossa ja sen jälkeen taulukon arvoille yritetään tehdä analyysiä matemaattisia menetelmiä ja teorioita käyttäen. Tällöin joudutaan arvioimaan funktion käyttäytymistä pisteissä, joista ei ole mittaustuloksia. Mikäli kyseessä on mittatuspisteiden välissä oleva piste, niin puhutaan interpoloinnista, muulloin suoritetaan  ekstrapolointia.

Taulukon muodostaminen

Ajatellaan, että mittaus on tuottanut lukupareja (xi,yi), i=1,...,5. Mathematica'ssa nämä parit kirjoitetaan listaksi

taulu = {{1, 3}, {2, 5}, {3, 7}, {4, 9}, {5, 12}}

TableForm[taulu]

Listan voi tulostaa taulukkomuodossa komennolla TableForm. Listan alkiot voi kertoa ja jakaa luvulla. Niistä voidaan vähentää ja niihin lisätä lukuja, mutta nämä laskutoimitukset kohdistuvat kaikkiin alkioihin. 
Esimerkiksi

taulu + 1

taulu*2

Jos pelkät yi:n arvot haluaa kertoa jollakin luvulla, niin silloin täytyy turvautua kahden listan kertolaskuun. Eli listan alkiot kerrotaan toisen listan vastinalkioilla. Esimerkkinä kerrotaan listassa taulu olevat yi:n arvot π:llä.

kerroin = Table[{1, Pi}, {i, 5}]

kerroin*taulu

Komennolla Table muodostetaan taulukko kerroin, jossa on viisi alkiota {1,π}. Yleisimmässä muodossa Table komennolla 
     Table[lauseke, {x, min, max, askel}]
muodostetaan taulukko, jossa x saa arvot parametrin min arvosta parametrin  max arvoon saakka, kun askeleen pituus on parametrin askel arvo. Lausekkeen lauseke arvo laskataan kullakin x:n arvolla ja nämä arvot muodostavat listan alkiot.
Taulukon alkioita voidaan käsitellä myös alkio kerrallaan. Kuten edellä mainittiin, yksittäiseen alkioon viitataan kaksinkertaisia hakasulkeita käyttäen. Esimerkiksi muodossa
	taulu[[i]]
Taulukon taulu y:n arvojen kertominen π:llä voisi tapahtua myös seuraavasti

maxalkio = Length[taulu]
Table[{taulu[[i, 1]], Pi*taulu[[i, 2]]}, {i, 
        maxalkio}]

Komento Length[taulu] kertoo kuinka monta alkiota taulukossa on.

Mathamaticassa voidaan datoja lukea myös tiedostosta ReadList-komennolla. 
Komennolla
	datat=ReadList["file", {Number,Number}]
luetaan tiedostosta file kahteen sarakkeeseen sijoitettuja lukuja taulukkoon datat. Lista {Number,Number} kertoo, missä muodossa taulukkon sarakkeet on talletettu tiedostoon.
Taulukoihin voidaan lisätä, poistaa ja muuttaa alkioita komennoilla 
Prepend[lista,alkio]       lisää alkio listan alkuun
Append[lista,alkio]        lisää alkio listan loppuun
Insert[lista,alkio,i]      lisää alkio listaan paikkaan i
Insert[lista,alkio,-i]     lisää alkio listaan paikkaan i lopusta lukien 
Delete[lista,{i,j,...]     poista alkiot i,j,... listasta
ReplacePart[lista,alkio,i] korvaa paikassa i oleva alkio

Interpolointi

Kun funktion arvo tunnetaan joukossa pisteitä, (yi=f(xi),i=1,...,n) ja halutaan laskea funktion arvo jossain muussa pisteessä x tarvitaan interpolointia.

Interpoloiva polynomi

Kun funktion arvo tunnetaan n:ssä pisteessä, voidaan sitä interpoloida n-1:n asteen polynomilla, joka kulkee kaikkien annettujen pisteiden kautta. Mathematicassa on käsky 
		InterpolatingPolynomial[datat,x], 
jolla tämä tehdään.
Testataan interpoloinnin tarkkuutta laskemalla sin(x) arvot pisteissä x=0,1,2,...,10 ja interpoloimalla saatua pistejoukkoa polynomilla.

datat = Table[{x, Sin[x]}, {x, 0, 10}] // N

TableForm[%]

poly = InterpolatingPolynomial[datat, x]

Expand[%]

poly /. x -> 1

Sin[1.5] - poly /. x -> 1.5

Sin[15] - poly /. x -> 15 // N

Polynomi on siis tarkka tunnetuissa pisteissä ja interpoloi pistejoukkoa melko tarkasti, mutta ekstrapolointia ei kannata käyttää tässä tapauksessa.

Palottain kuutiollinen sovitus

Käytännössä joudutaan interpolaatiota kuitenkin usein tekemään suurelle joukolle pisteitä, jolloin olisi käytettävä korkean asteen interpolaatiopolynomia. Laskennallisten ja laskentatarkkuuteen liittyvien ongelmien vuoksi on parempi käyttää palottaista alemman kertaluvun interpolaatiopolynomia. Palottain kuutiollinen sovitus saadaan komennolla Interpolation[data].

inter = Interpolation[datat]

Polynomin astetta voi muuttaa parametrilla {\tt InterpolationOrder->n}.

inter = Interpolation[datat, InterpolationOrder -> 2]

Tuloksena syntyi funktio, InterpolatingFunction[{{0,10.}},"<>"],  jonka määrittelyalue on {0, 10.}. Se sisältää riittävät tiedot interpoloivan funktion muodostamiseen kaikissa tarpeellisissa tilanteissa, kuten arvojen laskemisessa, piirtämisessä, derivoinnissa
jne.

inter[1.5]

Plot[inter[x], {x, 0, 10}]

Sin[1.5] - inter[1.5] // N

Huomaat, että Interpolation ja InterpolatingPolynomial eivät anna yhtä tarkkaa tulosta pisteessä 1.5.

Spline sovitus

Usein käytetty interpolointimenetelmä on sovittaa kolmannen asteen spline (cubic spline) datajoukkoon. Se määritellään polynomina
		s(x) = (Subscript[y, i])+(Subscript[b, i])(x-(Subscript[x, i]))+(Subscript[c, i])(x-(Subscript[x, i]))^2+(Subscript[d, i])(x-(Subscript[x, i]))^3
kun  (Subscript[x, i])<x<(Subscript[x, i + 1])
	- Derivaatat s'(x) ja s''(x) ovat jatkuvia kaikissa pisteissä (Subscript[x, i])
	- s((Subscript[x, i]))=(Subscript[y, i])
	- s((Subscript[x, i + 1]))=(Subscript[y, i + 1])
	- Funktio on mahdollisimman sileä, eli integraalilla
		\(TraditionalForm\`∫\_x1\%xn\(\((\(s'\)' \((x)\))\)\^2\) dx\) 
	  on minimi.
Spline interpolaatiota ei ole mathematican peruskäskyjen joukossa, mutta se löytyy paketista NumericalMath`SplineFit`. Lisäpaketeissa olevat käskyt saadaan käyttöön komennolla 
	 <<

<< NumericalMath`SplineFit`

spline = SplineFit[datat, Cubic]

spline[1.5]

Sovituksen tuloksena on funktio SplineFunction[Cubic, {0.,10.}, <>], jonka parametreina ovat pisteet numeroituna järjestyksessään. Funktion arvo esitetään (x,y) lukuparina.

Tehtäviä

24. Piirrä edellä muodostettua datajoukkoa datat interpoloivanpolynomin ja paloittain interpoloivan neliöllisen polynomin kuvaajat.
25. Muodosta erotusfunktio Sin[x]-IntepolatingPolynomial ja piirrä sen kuvaaja.

Kertaustehtäviä

4. Ratkaise yhtälö ln(x)=x-2.
5. Ratkaise lineaarinen yhtälöryhmä
	2 x+3 y-5 z=10
	13 x-z=1
	2y+z=3
6. Taulukoi lukuparit {x,ln x}, kun x=1,2,...,10. Muodosta paloittain interpoloiva toisen asteen polynomi ja piirrä sen kuvaaja.

Mittaustulosten käsittelyä

Pienimmän neliösumman sovitus

Erilaisia funktioita voidaan sovittaa mittaustuloksiin 
		Fit-komennolla. 
Esimerkkinä yksinkertainen pienimmän neliösumman suoran sovitus.

arvot = {{1, 1}, {2, 2}, {3, 4}, {4, 3}, {5, 4}, {6, 6}}

sovitus = Fit[arvot, {1, x}, x]

Piirretään nyt kuvat.

ListPlot[arvot]

Plot[sovitus, {x, 1, 6}]

Piirretään nyt kummatkin kuvat samaan koordinaatistoon.

Show[%, %%]

Myös muita funktioita voidaan sovittaa mittaustuloksiin.
 
Tilastollisen analyysin Fit-funktiolla tehdyn sovituksen hyvyydestä saa 
		Regress-funktiolla.
 Se löytyy Statistics`LinearRegression` paketista.
 
Esimerkiksi

<< Statistics`LinearRegression`

arvot = {{1, 1}, {2, 2}, {3, 2}, {4, 3}, {5, 5}, {6, 6}}

Regress[arvot, {1, x}, x]

antaa edellä käytettyihin datoihin arvot suoritetun lineaarisen sovituksen tilastollisen analyysin.

Tehtäviä

26. Sovita seuraavat mittausarvot toisen asteen polynomilla:
	(-2, 6),	(-1, 1),	(0, 0.5), (1, 2), (2, 4)
      Piirrä lisäksi mittausarvot ja sovitus samaan kuvaan.	
27. Tee äskeisen neliöllisen sovituksen tilastollinen analyysi Regress-funktiolla.
28. Kokeile sovittaa esimerkin suoraa toisen asteen polynomilla. Selitä mitä tapahtui.

Derivointi ja integrointi

Derivointi

Lausekkeen derivaatta saadaan komennolla 
          D[lauseke,muuttuja]
tai
          D[lauseke,{muuttuja,aste}]
Esimerkiksi

D[x^3  +  3 x^2  -  4 x  +  3, {x, 2}]

laskee lausekkeen (x^3) + 3(x^2) -4x + 3 toisen derivaatan x:n  suhteen.  Nopeus v on paikan x aikaderivaatta, v= \(TraditionalForm\`x\& . \)=dx/dt. Kiihtyvyys a puolestaan on nopeuden aikaderivaatta, eli paikan toinen aikaderivaatta, a= \(TraditionalForm\`\(x\& .. \)\)= \(TraditionalForm\`v\& . \)=dv/dt. Yläpisteellä merkitään aikaderivaattaa, \(TraditionalForm\`x\& . \)=dx/dt.

Tehtäviä

29. Derivoi x^2+4x.     
30. Derivoi x^4-3x-sin(x) kaksi kertaa.      
31. Kappaleen paikka saadaan lausekkeesta x=40-5t-5t^2. Mikä on
kappaleen kiihtyvyys hetkellä t=2?      
32. Piirrä edellisen tehtävän paikka ja nopeus ajan funktiona.    
33. Laske funktion f(x,y)=x^3 y^4
osittaisderivaatat
∂f / ∂x  ja  ∂f / ∂y   sekä  
(∂^2)f / ∂(x^2),   (∂^2)f / ∂(y^2),     (∂^2)f / ∂x∂y  ja   (∂^2)f / ∂y∂x .

Kokonaisdifferentiaali ja virheenarviointi

Laboratoriotöissä joudutaan usein tekemään  virheenarviointeja. Mathematica osaa laskea myös kokonaisdifferentiaalit näppärästi komennolla Dt.
       Dt[f]        laskee kokonaisdifferentiaalin df
       Dt[f,x]      laskee kokonaisderivaatan df/dx
       Dt[f,{x,n}]]  laskee astetta n olevan kokonaisderivaatan
       Tutkitaan lausekkeen 
       	f(a,b,c)=  \(TraditionalForm\`\(π\ \(a\^\(\(2\)\(\ \)\)\) b\^3\)\/\@\(\(c\)\(\ \)\)\)derivaattoja
ja kokonaisdifferentiaaleja.

Clear[a, b, c]
f  =  Pi a^2 b^3/Sqrt[c]

D[f, a]

D[f, {a, 2}]

Dt[f]

Dt[f, a]

Dt[f, {b, 2}]

Esimerkki

Halutaan tietää sylinterin tilavuus. Mitataan sylinterin halkaisija ja korkeus työntömitalla, kukin kolme kertaa. Tulokset olivat
Mittaus   Halkaisija, cm      Korkeus, cm
    1               12,00                       15,00
    2	        12,01                       14,95
    3                12,00                       15,05
Mittausten keskiarvot saadaan laskettua käyttämällä tilastofunktioiden pakkausta.

<< Statistics`DescriptiveStatistics`

data1  =  {12.00, 12.01, 12.00}

halkaisija  =  Mean[data1]

data2 = {15.00, 14.95, 15.05}

korkeus  =  Mean[data2]

d=halkaisija, cm         h=korkeus, cm
      12,0033                           15,00

Sylinterin tilavuus on
V = hπ ((((d/2)))^2)
\(ja\ mittausten\ antama\ virhearvio\ tilavuudelle\ saadaan\ laskemalla\)
\(halkaisijan\ ja\ korkeuden\ mittausvirheistä\ aiheutuneiden\ virheiden\)
\(itseisarvojen\ \(\(summa\)\(.\)\)


      \(\(|\)\(ΔV\)\(|\)\)\  = \ \(\(|\)\(\(∂\ 
                  V\/∂\ 
                  h\) Δh\)\(|\)\(+\ \(\(|\)\(\(∂\ 
                      V\/∂\ 
                      d\) Δd\)\(|\)\(
\)\(				\)\(
\)\(Tilavuuden\ kokonaisdifferentiaali\ lasketaan\ seuraavasti\)\)\)\)\)

Clear[h, d, v]
v[d_, h_] := h Pi (d/2)^2
diffV = Dt[v[d, h]]

virhe  =  Abs[diffV[[1]]]  + Abs[diffV[[2]]]

Sijoitetaan halkaisijan ja korkeuden arvoiksi mittausten keskiarvot ja differentiaalien paikalle niiden virheet.

N[virhe  /.  {Dt[d] -> 0.01,  Dt[h] -> 0.05,  d -> 12.0033,  
          h -> 15}]

Virhearvioksi saamme siten ±8.48

Tilavuus saadaan käyttämällä määriteltyä funktiota.

N[v[12.0033, 15]]

Nyt voimme kirjoittaa tuloksen V=1697±8 (cm^3).

Edellä differentiaalin numeroarvon laskemisessa käytettiin sijoitusoperaattoria  /. , jolla saatiin differentiaalilausekkeessa esiintyvät muuttujat korvattua numeroarvoilla.

Tehtäviä

34. Oletetaan, että eläintarhan apinoiden lukumäärä riippuu banaanien lukumäärästä ja yhden apinan painosta, kuten
m[banana_,paino_]=3*banana/paino
Laske apinoiden lukumäärä virherajoineen, kun banaaneita on 100 ± 5 kg ja yhden apinan paino on 7.3 ± 1.2 kg.

Integrointi

Mathematica osaa integroida melko hyvin jopa symbolisesti, mutta on muistettava että jos Mathematica ei osaa laskea esimerkiksi jotain integraalia symbolisesti, se ei tarkoita, etteikö integraalia voisi laskea käsin kohtuullisella vaivannäöllä. Komento
              Integrate[lauseke,x]
integroi annetun lausekkeen x:n suhteen. Komento
              Integrate[f,{x, a, b }]
laskee f:n  määrätyn integraalin  (∫_a^b) f(
            x) dx  ja komento 
              Integrate[f,{x, xmin, xmax},{y,ymin,ymax}]
laskee kaksinkertaisen integraalin. Numeerinen integrointi tapahtuu komennolla
              NIntegrate[f,{x, xmin, xmax}]
Numeerinen integrointi onnistuu paljon useammin kuin symbolinen integrointi, mutta se ei tietystikään ole aivan tarkkaa. Useimmissa tapauksissa Mathematica osaa kertoa, jos integroinnilla ei saavuteta haluttua tarkkuutta tai jos näyttää siltä, että integroitava funktio on liian hankala numeerisillekin menetelmille.

Tehtäviä

35. Integroi sin(x) välillä 0->2π. Kokeile myös numeerista integrointia.
36. Integroi (∫_a^b)(∫_d^c) 
                (1/(x y)) dx dy.
37. Integroi ∫  (4 x^4 + 4 x^3 -9 x ^2 -x +2 ) dx  ja tarkista saamasi tulos derivoimalla se.
38. Ylinopeutta ajava auto ajaa nopeudella v=15 m/s.  Poliisi ryhtyy ottamaan ajajaa kiinni. Poliisiauto lähtee levosta ja kiihdyttää a=10 m/s^2.  Milloin poliisi saa kaaharin kiinni?  Piirrä kummankin auton paikka ajan funktiona.
39. Kappaleen kiihtyvyys a(t) ajan funktiona muuttuu seuraavasti
	a(t) = arctan(t),
kun aika ilmoitetaan sekunteina. Laske, mikä on kappaleen keskimääräinen nopeus 12 ensimmäisen sekunnin aikana, kun liikkeelle lähdetään levosta, ja mikä on sen nopeus 20 sekunnin kuluttua lähdöstä. Piirrä kiihtyvyyden ja nopeuden kuvaajat samaan kuvaan.

Raja-arvot  ja sarjat

Raja-arvot

Raja-arvojen laskeminen Mathematicalla tapahtuu komennolla 
          Limit[f[x],x->x0]
Tässä lasketaan raja-arvo (lim_(x -> x0))f(x).
Esimerkiksi funktion sin(x)/x raja-arvo nollassa:

f = Sin[x]/x

Limit[f, x -> 0]

Taylorin ja Laurentin sarjat

Mathematica osaa kehittää funktioita potenssisarjaksi. Taylorin ja Laurentin sarjakehitelmät saadaan Series-komennolla.		
Esimerkiksi sinin sarjakehitelmä muuttujan x suhteen pisteen x=0 ympäristössä potenssiin yhdeksän saakka  saadaan komennolla:

Series[Sin[x], {x, 0, 9}]

Normal[%]

Viimeinen O[x]-termi kuvaa yhdeksättä potenssia korkeampia termejä. Ennen kuin Mathematican tekemää sarjakehitelmää voidaan käyttää edelleen (esimerkiksi kuvaajan piirtoon), pitää se muuttaa tavalliseksi Mathematican lausekkeeksi Normal-komennolla.

Laurentin sarjassa on termejä, jotka ovat äärettömiä siinä pisteessä, jonka ympäristössä kehitelmä on tehty. Esimerkiksi funktion ((e^x)/(1 - (x^2)))Laurentin sarjakehitelmä pisteen x=1 ympäristössä on

Series[E^x/(1 - x^2), {x, 1, 3}]

Tehtäviä

40. Kehitä kosini sarjaksi yhdeksänteen asteeseen asti. Vertaa sitä esimerkin sinin sarjakehitelmään.
41. Laske yhteen sinin ja kosinin sarjakehitelmien neliöt. Mitä pitäisi tulla tulokseksi ja mitä saat?
42. Piirrä sinin kolmannen, viiden, seitsemännen ja yhdeksännen asteen sarjakehitelmien kuvaajat ja vertaa niitä tarkkaan kuvaajaan.
43. Kehitä sini sarjakehitelmäksi pisteen x=1 ympäristössä kuudennen asteen termejä myöten.Esitä tulos polynomina käyttäen Expand-komentoa.

44. Tutki myös jonkin muun funktion sarjakehitelmien tarkkuutta muunnellen astelukua.
45. Kokeile kehittää jokin määrittelemätön funktio f[x] sarjakehitelmäksi.Mitä saat tulokseksi?
46. Esimerkkinä singulaarisista eli ns.Laurentin sarjakehitelmistä kehitä funktio (e^x)/(x-1) sarjaksi pisteen x=1 ympäristössä.Mathematican komento on sama kuin Taylorin sarjakehitelmän tapauksessa.
47. Muunna kehitelmä Normal-funktion avulla normaaliin lausekkeen muotoon kuvan piirtämistä
varten ja piirrä kuva.
48. Osoita,että suhteellisuusteorian mukaisesta energian lausekkeesta E=m(c^2)saadaan ei-relativistinen liike-energian lauseke
E=(Subscript[m, 0])(c^2)+1/2(Subscript[m, 0])(v^2),
kun liikemassan m ja lepomassan (Subscript[m, 0]) välille voidaan johtaa relaatio m=(Subscript[m, 0])/sqrt(1-v^2/c^2).

Differentiaaliyhtälöitä Mathematica'lla

Differentiaaliyhtälöt

Differentiaaliyhtälöt ovat yhtälöitä, joissa esiintyy tuntematon funktio, esimerkiksi y(x), ja sen derivaattoja.
Yksinkertainen esimerkki differentiaaliyhtälöstä on

		((dy(x))/dx)= 3y(x)

Differentiaaliyhtälön ratkaisuksi sanotaan analyyttistä lauseketta tai  muuta esitystä y(x):lle. Yhtälön ratkaisu on y(x)=(e^(3  x)), sillä

		((d (e^(3  x)))/dx)=3 (e^(3  x))= 3y(x) 

Differentiaaliyhtälön kertaluku on korkein yhtälössä esiintyvä derivaatan kertaluku. Yo. yhtälö  oli ensimmäisen kertaluvun differentiaaliyhtälö. Esimerkki toisen kertaluvun differentiaaliyhtälöstä on

	(((d^2) u(x)  )/((d^2) x))+16 u(x) = 0

sillä siinä esiintyy toista derivaattaa. Tämän yhtälön ratkaisu on mm. u(x) = cos(4x).

Reunaehdot ja reuna-arvotehtävät

Tarkastele seuraavaa differentiaaliyhtälöä, joka kuvaa populaation määrää x ajanhetkellä t.
	\(TraditionalForm\`x\& . \)= (dx/dt)= x
Tämän yhtälön ratkaisu on x(t)=k(e^t). Koska ratkaisussa esiintyy mielivaltainen vakio k, kutsutaan tällaista ratkaisua yleiseksi ratkaisuksi. Jos toisaalta tiedämme, että ajanhetkellä t=0 x(0)=10, voimme sijoittaa nämä tiedot yleiseen ratkaisuun ja saamme k:lle arvon k=10. Olemme juuri ratkaisseet reuna-arvotehtävän.

Ratkaisut

Jos differentiaaliyhtälö ratkaistaan ilman reunaehtoja, saadaan ratkaisuun yhtälön kertaluvun osoittama määrä tuntemattomia vakioita. Ratkaisu on siis käyräparvi, jonka parametrien määrä on yhtälön kertaluku.

Lineaarisia differentiaaliyhtälötyyppejä

Lineaarisuus tarkoittaa, että yhtälössä ei esiinny tuntemattoman funktion y(x) tai sen derivaattojen potensseja. Esittelemme muutamia differentiaaliyhtälötyyppejä, joihin tulet törmäämään fysiikassa.

Ensimmäisen kertaluvun vakiokertoimiset differentiaaliyhtälöt

Ensimmäisen kertaluvun vakiokertoimiset differentiaaliyhtälöt ovat  muotoa
		y'(x)+ a y(x)=c,
missä a ja c ovat vakioita. Jos c=0 sanotaan yhtälöä homogeeniseksi.

Ensimmäisen kertaluvun funktiokertoimiset differentiaaliyhtälöt

Muotoa
		y'(x)+p(x)y(x)=q(x)
olevat differentiaaliyhtälöt ovat ensimmäisen kertaluvun funktiokertoimisia differentiaaliyhtälöitä.  Jos q(x)=0 kaikilla x:n arvoilla on yhtälö homogeeninen, jos taas sekä q(x) että p(x) ovat nollasta eroavia, sanotaan yhtälöä täydelliseksi.

Toisen kertaluvun vakiokertoimiset
differentiaaliyhtälöt

Toisen kertaluvun differentiaaliyhtälöissä esiintyy jo tuntemattoman funktion y(x) toista derivaattaa y''(x). Vakiokertoiminen, täydellinen yhtälö on muotoa
		y''(x)+a y'(x) + b y(x)=c.
Vastaavasti kuin ensimmäisen kertaluvun tapauksessa sanotaan yhtälöä homogeeniseksi, jos c=0.

Toisen kertaluvun funktiokertoimiset
differentiaaliyhtälöt

Yleinen muoto toisen kertaluvun funktiokertoimiselle differentiaaliyhtälölle on
	y''(x)+p(x) y'(x) + r(x) y(x)=q(x).

Differentiaaliyhtälöiden ratkaiseminen Mathematicalla

Mathematica on näppärä apuneuvo differentiaaliyhtälöiden ratkaisuun.
Peruskomento on DSolve.
	DSolve[yhtälö,y[x],x]
Ratkaisemme esimerkiksi yhtälön

DSolve[y'[x] == 3 y[x], y[x], x]

Tulosta pitää osata hieman tulkita, mutta kun hoksaa että siinä esiintyvä C[1] on vakio, jota olimme aiemmassa esimerkissä merkinneet k:lla, lauseke näyttää jo tutulta.  Huomaa, että jos erehdyt käyttämään  =-merkkiä ==-merkin asemesta, määrittelee Mathematica  =-merkin vasemmalla puolella olevan asian oikealla puolella olevan asian kanssa yhtäsuureksi. Jos näin tapahtuu esimerkiksi ratkottaessa differentiaaliyhtälöitä, pitää nämä määrittelyt purkaa Clear-komennolla, muuten voi syntyä ihmeellisiä virheitä.

Numeerinen ratkaisu

Jos differentiaaliyhtälö on hankala, ei Mathematica välttämättä osaa ratkaista sitä symbolisesti. Voit kuitenkin käyttää komentoa 
NDSolve, joka ratkaisee differentiaaliyhtälön
numeerisesti.  NDSolved:lle pitää antaa myös vapaan muuttujan väli, jolle ratkaisu etsitään. Tulos on voimassa vain tällä välillä. Lisäksi NDSolve vaatii riittävästi reunaehtoja kyetäkseen ratkaisemaan kaikki vakiot, muuten se ei pysty laskemaan numeerista ratkaisua.
Esimerkiksi edellisen esimerkin yhtälön tapauksessa NDSolve ei antaisi tulosta, mutta jos sovimme reunaehdoksi vaikkapa y(0)=1 saamme,

ratkaisu = 
      NDSolve[{y'[x] == 3 y[x], y[0] == 1}, y[x], {x, 0, 10}]

Nyt voimme esimerkiksi piirtää funktion y(x) kuvaajan.

Plot[y[x]  /.  ratkaisu, {x, 0, 3}]

Kuvaaja voidaan piirtää vain välillä [0,10]. Huomaa sijoitusoperaattorin /. käyttö!

Tehtäviä

49. Ratkaise 1:n kertaluvun lineaarinen vakiokertoiminen differentiaaliyhtälö y'(x)-3y(x)=4. 
50. Ratkaise 2:n kertaluvun lineaarinen vakiokertoiminen differentiaaliyhtälö 4 y''(x)-3y'(x)+5y(x)=2. Määrää mieleisesi reunaehdot ja piirrä ratkaisun reaaliosan kuvaaja.

Kertaustehtäviä

7. Integroi cos(x) välillä x∈(-10,10) symbolisesti ja numeerisesti.
8. Kirjoita tan(x):n sarjakehitelmä 5:n ympäristössä. Ota kehitelmään 5 ensimmäistä termiä.
9. Ratkaise differentiaaliyhtälö 12 x''(t)-9 x'(t)+16 x(t)-9=0 alkuarvoilla x(0)=2 ja x'(0)=0. Piirrä ratkaisun kuvaaja, kun t saa arvot 0->10.
10. Mieti, miten saisit näppärästi määriteltyä edellisen tehtävän ratkaisua vastaavan funktion.

Radioaktiivinen hajoaminen

Radioaktiivista hajoamista voidaan kuvata yhtälöllä
(dN/dt)= -λ N
missä  λ on hajoamisvakio.

Hetkellä t =0 meillä on (Subscript[N, 0])  ydintä. Hetkellä t jäljellä olevien ydinten määrä saadaan selville ratkaisemalla yo. yhtälö 
alkuehdolla N(t=0)=(Subscript[N, 0]).

ratk =  
      DSolve[{n'[t] == -lambda n[t], n[0] == n0}, n[t], 
        t]

joka voidaan kirjoittaa tutumpaan muotoon
	N(t)=(Subscript[N, 0])(e^(-λt)).
Lasketaan vielä puoliintumisaika. Puoliintumisaika on aika, jonka kuluttua puolet ytimistä on hajonnut.

Solve[n[t] == 0.5 n0  /.  ratk, t]

Voimme myös havainnollistaa ratkaisua piirtämällä siitä kuvan. Ensin annamme λ:lle ja (Subscript[N, 0]):lle sopivat numeroarvot.

n0 = 1000

lambda = 0.03

Plot[n[t]  /.  ratk, {t, 0, 100}]

Huomaa, miten ratkaisuna saamaamme funktioon n[t] viitataan sekä Solve- että  Plot-komennoissa.

Ylöspäin heitetty kappale

Haluamme saada matemaattisen kuvauksen suoraan ylöspäin heitetylle pesäpallolle. Tiedämme, että pallon kiihtyvyys on gravitaatiokiihtyvyys. Toisaalta kiihtyvyys on paikan toinen aikaderivaatta, 
	a=g=(((d^2) y)/(dt^2))= \(TraditionalForm\`\(y\& .. \)\).
Lisäksi voimme asettaa koordinaatiston niin, että pallo on lähtöhetkellä origossa.  Matemaattisesti tämä voidaan esittää yhtälöparilla
	a=\(TraditionalForm\`\(y\& .. \)\)
	y(t=0) = 0
Katsotaan mitä Mathematica saa näistä aikaiseksi.

DSolve[{a == y''[t], y[0] == 0}, y[t], t]

Huomaa, että olimme merkinneet \(TraditionalForm\`\(y\& .. \)\):n t:n funktioksi, vaikka kyseessä oikeastaan olikin vakio a. Huomaa myös, että tuloksessa on ainoastaan yksi tuntematon vakio, C[2], vaikka kyseessä oli toisen kertaluvun differentiaaliyhtälö. Tämä johtuu siitä, että annoimme yhtälön ratkaisun yhteydessä reunaehdon, y(0)=0. Pääsisimme eroon myös C[2]:sta, jos meillä olisi vielä yksi lisäehto, esimerkiksi tieto siitä, että pesäpallon alkunopeus oli 10 m/s.

Tehtäviä

51. Eliminoi esimerkin tehtävästä vakio C[2] käyttämällä oletusta, jonka mukaan pesäpallon alkunopeus oli 10 m/s ylöspäin.
52. Mikä on pesäpallon saavuttama maksimikorkeus ja missä ajassa? Milloin pesäpallo palaa heittokorkeudelle?

Putoamisliike vedessä

Tutkitaan nyt hieman monimutkaistettua tilannetta, jossa tiputamme tiiliskiven veneestä. Tiiliskiveen vaikuttaa maan vetovoima F_G=mg, mutta myös veden aiheuttama liikevastus F_D=-k\(TraditionalForm\`y\& . \). Muotoillaan asia matemaattisesti.
		F = mg-k \(TraditionalForm\`y\& . \) = m \(TraditionalForm\`\(\(y\& .. \)\(\ \)\)\)= ma
Mathematica sanoo

Clear[y, m, g, t]
DSolve[m g  -  k y'[t] == m y''[t], y[t], 
      t]

Kun otetaan huomioon järkevät reunaehdot, saadaan tulokseksi

DSolve[{m g  -  k y'[t] == m y''[t], 
        y[0] == 0, y'[0] == 0}, y[t], t]

Tehtäviä

53. Päättele veneestä pudotetun tiiliskiven tapauksessa järkevät reunaehdot, ratkaise yhtälö Mathematicalla, anna vakioille numeroarvot ja piirrä ratkaisusta kuva. Mitä huomaat? Miten tutkisit asiaa edelleen? Tutki myös mitä tapahtuu, jos muutat tiiliskiven massaa,ilman että muutat vakiota k.
54. Laske tiiliskiviesimerkki numeerisesti.Vertaa tuloksien yhtäpitävyytta piirtämällä kummastakin kuvaajat päällekäin, tai vaikka laskemalla niiden arvot muutamassa pisteessä.

Raketin lento

Tarkastellaan rakettia, joka polttaa massaansa tuottaakseen työntövoimaa. Rakettiin vaikuttaa vain maan vetovoima -m(t)g, joka oletetaan vakioksi, ja raketin moottorin aiheuttama työntövoima (Subscript[F, moot]), joka myös oletetaan vakioksi. 
Liikeyhtälön muodostaa yhtälöpari
	(Subscript[F, moot])-m(t)g=m(t) \(TraditionalForm\`\(x\& .. \)\)(t)
	m(t)=(Subscript[m, 0])-a t
niin kauan kuin poltettavaa massaa riittää. Kun polttoaine on kulutettu muuttuu m(t):n lauseke. Tutkitaan nyt kuitenkin rakettimatkan alkuosan t<10sek ratkaisua seuraavalla ohjelmanpätkällä.

Clear[m]
g = 9.81
m0 = 100
a = 5
f = 3000
m[t_] = m0 - a t
\(\(\(ratkaisu = 
      NDSolve[{f\  - \ m[t]\ g == m[t]\ \(\(x'\)'\)[t], 
          x[0] == 0, \(x'\)[0] == 0}, x[t], {t, 0, 10}]\)\(
\)
    \)\)

Koska tulos ei ole kovin selväkielinen, tutkimme asiaa piirtämällä ratkaisusta, eli raketin korkeudesta ajan funktiona, kuvan.

Plot[x[t]  /.  ratkaisu, {t, 0, 10}]

Tehtäviä

55. Laske rakettiesimerkki ja piirrä sen ratkaisu,mutta lisää ilmanvastuksen aiheuttama voima,  (Subscript[F, ilma])=-kv^2=-k \(TraditionalForm\`x\& . \)(t)^2. Käytä k:lle vaikka arvoa 0.3. Vieläkö ratkaistava yhtälö on lineaarinen? Vertaa raketin lentorataa esimerkin lentorataan.

Mekaniikkaa

Kitkaton värähtelijä

Tutkitaan kitkattoman jousen päähän kiinnitetyn kappaleen käyttäytymistä. Jousen jousivakio on k ja  kappaleen massa m. Systeemi oletetaan kitkattomaksi. Kappaletta poikkeutetaan matkan a verran tasapainoasemastaan ja päästetään vapaaksi. Haluamme tietää, mikä on punnuksen poikkeama x tasapainoasemastaan  kullakin ajanhetkellä. Kirjoitamme näkyviin liikeyhtälön:
	F= ma=-kx
eli	
	 m \(TraditionalForm\`\(x\& .. \)\) = -kx
Lisäksi tiedämme, että punnus oli hetkellä t = 0 poikkeutettuna matkan a tasapainoasemastaan, x(0) = a ja sen nopeus oli nolla, 		
	(OverDot[x])(0)=v(0)=0. 
Voimme käyttää näitä tietoja reunaehtoina.
Syötämme tietomme Mathematicaan:

Clear[m, x, k, a]
DSolve[{m x''[t] == -k x[t], x[0] == a, x'[0] == 0}, 
      x[t], t]

Tehtäviä

56. Laske esimerkin jouseen kiinnitetyn kappaleen hetkellinen nopeus.
57. Kirjoita tiedostoon Mathematica-ohjelma, joka antaa a:lle, k:lle ja m:lle järkevät arvot,ratkaisee yhtälön kuten edellä ja piirtää kuvaajan kappaleen poikkeamasta ajan funktiona.
58. Lisää ohjelmaan pätkä, joka laskee kappaleen hetkellisen nopeuden ja piirtää sen samaan kuvaan poikkeaman kanssa.Tulkitse kuvaa. (Huomaa, ettet voi viitata tuloksiin %n-menetelmällä ohjelman sisällä!)

Vaimennettu värähtely

(Katso Benson, University Physics s.309, 15.5)
Todellisuudessa kuvan  kaltaisessa systeemissä esiintyy erilaisia kitkavoimia. Kun jouseen kiinnitetyn kappaleen nopeus on pieni ja systeemi upotetaan nesteeseen, on kitkavoima verrannollinen nopeuteen,
	(Subscript[F, kitka])=-γ v.
Suuretta γ sanotaan vaimennusvakioksi.  Kirjoitamme uudelleen liikeyhtälön ja otamme huomioon myös kitkan
		m \(TraditionalForm\`\(x\& .. \)\)+γ\(TraditionalForm\`x\& . \) +kx = 0

Tehtäviä

59. Ratkaise jousen paikka ajan funktiona.
60. Kuten kitkattomassa tapauksessa, piirrä paikan ja nopeuden aikakehitys samaan kuvaajaan. Mitä eroja havaitset?

Sähköoppia

LC-piirit

Tarkastellaan suljettua virtapiiriä, joka koostuu kondensaattorista, kapasitanssi C  ja kelasta, induktanssi L. Kirchhoffin lain mukaan kierrettäessä virtapiiri, tulee jännitteiden (potentiaalierojen) summan olla nolla, (Subscript[V, C])+(Subscript[V, L])=0. Sähköopista tiedetään että virta on sähkövarauksen aikaderivaatta, I=dQ/dt. Keloille pätee (Subscript[V, L])=L dI/dt, eli (Subscript[V, L])=L (d^2)Q/(dt^2). Kondensaattorin napojen välillä vallitsee jännite (Subscript[V, C])=Q/C

		(Subscript[V, L])=L(d^2)Q/(dt^2)
		(Subscript[V, C])= Q/C
		(Subscript[V, C])+(Subscript[V, L]) = 0
		Q/C+L(d^2)Q/(dt^2)=0

Yllättäen edessämme on toisen kertaluvun differentiaaliyhtälö. Ei kun ratkaisemaan! Mathematica sanoo

Clear[c, l, q]
DSolve[q[t]/c + l q''[t] == 0, q[t], t]

Saimme harmonisen oskillaattorin lausekkeen.

Tehtäviä

61. Eliminoi vakiot asettamalla varaukselle q ja virralle q' sopivat alkuehdot, ja tutki saamaasi lopputulosta kuvin.Vaihtele komponenttien arvoja ja tutki systeemin käyttäytymistä.

LCR-piirit

LRC-piiri on huomattavasti lähempänä todellista tilannetta. Siinä on kytketty sarjaan kondensaattorin ja kelan kanssa myös vastus, joka aiheuttaa sähköenergian muuttumisen lämpöenergiaksi ja häviämisen systeemistä.  Kirchhoffin säännön perusteella voimme kirjoittaa yhtälön
	
	Q/C+IR+LdI/dt=0
	
Käyttämällä LC-piirien yhteydessä tekemäämme havaintoa I=dQ/dt voimme kirjoittaa saman muodossa
 	
 	L\(TraditionalForm\`\(Q\& .. \)\)+R\(TraditionalForm\`Q\& . \)+Q/C=0
 	
 Verrattuna edelliseen tehtävään olemme saaneet yhden termin lisää. Mathematica osaa ratkaista tämän yhtälön analyyttisesti, mutta tulos on melko sekavan näköinen

DSolve[{l q''[t]  +  r q'[t]  +  q[t]/c  ==  0}, 
      q[t], t]

Simplify[%]

Tehtäviä

62. Eliminoi vakiot kuten LC-piirien tapauksessa, ja tutki saamaasi lopputulosta kuvin. Vaihtele komponenttien arvoja ja  tutki systeemin käyttäytymistä.

Sähköopin ja mekaniikan suureiden analogiaa

Mekaniikka  x     v      m      m(v^2)/2       k        k(x^2)/2      F

Sähköoppi   Q    I        L          L((I )^2)/2   1/C   (Q^2)/2C     V

Populaatioiden kilpailu

Epälineaarisia differentiaaliyhtälöitä

Edellä on tutkittu lineaarisia differentiaaliyhtälöitä, joilla on paljon sovellutuksia eri fysiikan aloilla. Seuraavassa tarkastellaan eräitä epälineaarisia differentiaaliyhtälötyyppejä, joilla voidaan kuvata toisiinsa kytkettyjen systeemien kehittymistä. Otamme tutkimuksen kohteeksi eläin- ja kasvimaailman eliöstöjen lisääntymisen ongelmat. 
Perusyhtälö lukumäärän N(t) kehittymiselle ajan t funktiona saadaan, jos tunnetaan lisääntymisnopeus ϵ,

		((dN(t))/dt)= ϵ N(t).
Yhtälö on muodoltaan sama kuin radioaktiivista hajoamista kuvaava yhtälö.  Se on lineaarinen yhtälö ja tuloksena on eksponentiaalinen kasvu ilman rajoja. Biologisen systeemin käyttäytyminen tällä tavoin on mahdollista vain tiettyyn  rajaan saakka, sillä kasvua rajoittavat tilan, ruuan, valon jne puute, eli yksilöt tulevat tietoisiksi ympäristöstään ja lajitovereidensa olemassaolosta. Tämä näkyy lukumäärän kehittymistä kuvaavassa differentiaaliyhtälössä epälineaarisena terminä, jonka edessä kertoimena on erilaisia puutteita kuvaava vakio β.

		((dN(t))/dt) = ϵ N(t) - β (N^2)(t) .

e  =  4
b   =  2
ra = 
      DSolve[{n'[t] == e n[t] - b n[t]^2, n[0] == 1}, n[t], 
        t]

Plot[n[t] /. ra, {t, 0, 10}, PlotRange -> All]

Toisenlainen epälineaarinen kytkentä aiheutuu toisten lajien olemassaolosta. Peto-saalis suhteessa saaliin määrä lisääntyy eksponentiaalisesti, ellei petoja ole olemassa. Toisaalta petojen lisääntyminen riippuu täysin saaliin lukumäärästä ja ne kuolevat pois eksponentiaalisesti, mikäli saalis hävitetään. Tällaista tilannetta kuvaamaan tarvitaan kaksi kytkettyä differentiaaliyhtälöä. Olkoon esimerkkinä kettujen ja jänisten lukumäärien kehittyminen,

	((d(Subscript[N, j]))/dt)= [ ϵ\(TraditionalForm\`\_j\) -  γ\(TraditionalForm\`\_j\) (Subscript[N, k])(t) ]  (Subscript[N, j])(t)
  
  	((d(Subscript[N, k]))/dt)= [- ϵ\(TraditionalForm\`\_k\) -  γ\(TraditionalForm\`\_k\) (Subscript[N, j])(t) ]  (Subscript[N, k])
  
Yhtälöissä (Subscript[N, j])(t) kuvaa jänisten ja (Subscript[N, k])(t) kettujen lukumäärää. Vakio ϵ \(TraditionalForm\`\_j\) kertoo jänispopulaation kasvunopeuden ilman kettuja. Se riippuu ruuasta, ilmastosta jne. Vakio γ\(TraditionalForm\`\_j\)  kuvaa sitä, kuinka paljon jäniksiä joutuu kettujen suihin, ja γ\(TraditionalForm\`\_k\)puolestaan sitä, miten kettujen määrä lisääntyy tämän saalistuksen seurauksena. Ketut eivät elä ilman jäniksiä, joten ϵ \(TraditionalForm\`\_\(\(k\)\(\ \ \)\)\)  antaa nopeuden niiden sukupuuttoon kuolemiselle. Kokeilemalla erilaisia parametrien arvoja yhtälön ratkaisemiselle voit todeta, että luonnon tasapaino on hyvin herkkä ja sen järkkyminen johtaa aivan uudenlaiseen ratkaisuun.	
 Mathematica'ssa nämä yhtälöt saavat muodon,

e1 = 10
g1 = 0.1
e2 = 0.6
g2 = 0.1

\(rat = 
      NDSolve[{\(j'\)[t] == \((e1 - g1\ k[t])\) j[t], \ \(k'\)[
              t] == \((\(-e2\) + g2\ j[t])\) k[t], 
\ \ \ \ \ \ \ \ j[0] == 
            1, k[0] == 100}, {j[t], k[t]}, {t, 0, 10}]\)

Plot[{j[t] /. rat, k[t] /. rat}, {t, 0, 10}, PlotRange -> All, 
      PlotStyle -> {Hue[0], Hue[0.5]}]

\(\(\(ParametricPlot[{j[t] /. rat[\([1]\)], k[t] /. rat[\([1]\)]}, {t, 0, 
        1}, PlotRange -> All]\)\(
\)\)\)

Yllä esitetyt mallit ovat tietysti hyvin idealisoituja ja luonnossa monet tekijät yhdessä vaikuttavat lukumäärien kehittymiseen. Esimerkkinä voisi olla mäntyjen ja kuusien kilpailu elintilasta. Sitä voidaan kuvata differentialliyhtälöryhmällä, jossa molemmat edellä esitetyt epälineaariset termit on otettu mukaan. 

	((d(Subscript[N, m]))/dt)= [ ϵ\(TraditionalForm\`\_m\) - β\(TraditionalForm\`\_m\) (Subscript[N, m])(t) - γ\(TraditionalForm\`\_m\) (Subscript[N, k])(t) ]  (Subscript[N, m])(t)
  
  	((d(Subscript[N, k]))/dt)= [ϵ\(TraditionalForm\`\_k\)  -  β\(TraditionalForm\`\_k\) (Subscript[N, k])(t)   -  γ\(TraditionalForm\`\_k\) (Subscript[N, m])(t)] (Subscript[N, k])(t)

Yhtälöissä (Subscript[N, m])(t) kuvaa mäntyjen ja (Subscript[N, k])(t) kuusien lukumäärää. Vakiot ϵ\(TraditionalForm\`\_m\) ja  ϵ\(TraditionalForm\`\_k\) kertovat  puiden lukumäärän lisääntymisnopeuden. Lukumäärän kasvaessa alkavat toiset puut varjostaa ja estää toisten kasvua, joten lukumäärä ei pääse kasvamaan rajatta. Vakiot β\(TraditionalForm\`\_m\) ja β\(TraditionalForm\`\_k\)  kuvaavat saman lajin sekä γ\(TraditionalForm\`\_m\)  ja γ\(TraditionalForm\`\_k\)  eri lajin aiheuttamaa varjostusta.
Kuva lukumäärän kehittymisestä ajan funktiona saadaan helposti piirtämällä Plot-funktiolla yhtälöiden ratkaisut. Toinen mielenkiintoinen tapa yhtälöryhmän ratkaisujen tutkimisessa on esittää ratkaisut toistensa funktiona, eli kuinka kettujen lukumäärä muuttuu jänisten lukumäärän funktiona. Tähän tarkoitukseen on käytettävissä ParametricPlot-funktio.  Se on muotoa 

	ParametricPlot[{(Subscript["\!\(\*StyleBox[\"\\\"f\\\"\", \"TI\"]\)", x]), (Subscript["\!\(\*StyleBox[\"\\\"f\\\"\", \"TI\"]\)", y])}, {t, tmin, tmax}] 

Tässä  piirretään funktion f(t) x- ja y-komponentit  (Subscript["\!\(\*StyleBox[\"\\\"f\\\"\", \"TI\"]\)", x]) (t) ja (Subscript["\!\(\*StyleBox[\"\\\"f\\\"\", \"TI\"]\)", y])

e1 = 1
g1 = 0.1
b1  =  0.1
e2 = 0.6
g2 = 0.1
b2  =  0.1

dy1  :=  m'[t] == (e1 - b1 m[t] - g1 k[t]) m[t]
dy2  :=  k'[t] == (e2 - b2 k[t] - g2 m[t]) k[t]
rat = 
      NDSolve[{dy1,  dy2, m[0] == 1, k[0] == 100}, {m[t], k[t]}, {t, 0, 
          100}]

Plot[{m[t] /. rat, k[t] /. rat}, {t, 0, 100}]

\(\(\(ParametricPlot[{m[t] /. rat[\([1]\)], k[t] /. rat[\([1]\)]}, {t, 0, 
        100}, PlotRange -> All]\)\(
\)\)\)

Tehtäviä

63. Tutki yhtälön dN(t)/dt= ϵ N(t) - β(N^2)(t) käyttäytymistä eri parametrien arvoilla ja antamalla lukumäärälle eri lähtöarvoja ajan hetkellä t=0. Piirrä lukumäärän kuvaaja ajan funktiona.
64. Tutki sekä jänisten ja kettujen että kuusien ja mäntyjen elämää kuvaavien yhtälöiden käyttäytymistä eri parametrien arvoilla ja eri lähtöarvoja käyttäen. Piirrä kettujen lukumäärä jänisten lukumäärän funktiona ja kuusien lukumäärä mäntyjen lukumäärän funktiona. Etsi kuvaajien avulla stabiilit ratkaisut, joissa lajien lukumäärät pysyvät ajan funktiona vakioina.

Ohjelmointi Mathematica'lla

Ohjelma ja sen suoritus

Mathematica'n kieli soveltuu myös ohjelmointiin. Se tuntee yleisimmät toistorakenteet ja ehtolauseet. Mathematica-ohjelmat kirjoitetaan tavalliseen tekstitiedostoon, jonka tunnisteeksi on sovittu .m . Esimerkiksi ohjelma.m on hyvä tiedostonimi. 
Mathematican kieli on tulkittu kieli, toisin kuin esimerkiksi C-kieli Atk I-kurssilla. Tämä tarkoittaa sitä, että Mathematica lukee ohjelmaa ja käy sitä suoraan läpi käsky käskyltä, kääntämättä sitä ensin mihinkään muuhun - kuten esimerkiksi assembly - muotoon. 
Kun ohjelma on kirjoitettu tiedostoonsa, käsketään Mathematicaa suorittamaan se. Mathematica lukee ohjelman käsky käskyltä ja yrittää suorittaa sen. Kaikkien grafiikkakäskyjen tulosteet avautuvat myös grafiikkanäytölle.
Ennen ohjelman ajoa on syytä varmistaa, että $Path-määrittelyissä on mukana hakemisto, jossa ohjelmatiedosto on. Hakemisto liitetään määrittelyihin komennolla SetDirectory["hakemisto"].
Ohjelma luetaan Mathematicaan << - käskyllä, esimerkiksi

Clear[a, b, c, f]
<<  ohjelma.m

Ohjelma oli

Clear[f, x]
a = 1.5
b = -2
c = 3
f[x_] := a x^2  +  b x  +  c
\(Solve[f[x] == 0, x]
\)

Tehtäviä

65. Tutki esimerkin ohjelmaa ohjelma.m. Käy se läpi käsky käskyltä, ja selitä, mitä se tekee.
66. Kirjoita itse tiedostoon ohjelma, jossa määrittelet funktion f(x)=x^3-4x+12 ja piirrät sen kuvaajan.

Ehtolauseet

Valinta eri käskyjen suorittamisen välillä tehdään If-komennolla
            If[ehto
          	, käskyt1
               , käskyt2]

Siinä annetaan ehto, jonka perusteella valinta suoritetaan. Ehto voi saada kaksi arvoa tosi ja epätosi. Ehdon jälkeen seuraa kaksi käskysarjaa pilkuilla erotettuna. Sarja käskyt1 suoritetaan, jos ehto on tosi ja  käskyt2 jos ehto on epätosi. On myös mahdollista antaa kolme lauseketta, joista viimeinen suoritetaan, jos ehto ei ole sen paremmin tosi kuin epätosikaan.

If[1 < 2, Print["Tosi"], Print["Vaarin"]]

Ehtorakenne, joka usein tulee vastaan matemaattisia funktioita käsiteltäessä, liittyy funktion määrittelyalueen rajaamiseen. Funtioiden ei tarvitse olla edes jatkuvia, vaan ne voidaan määritellä paloittain. 
	funktio[x_]:=lauseke1/;määrittelyalue1
        	funktio[x_]:=lauseke2/;määrittelyalue2
Merkintää  /;  käytetään siis määrittelyalueen rajaamiseen.

Vertailuoperaattorit

Ehtojen muodostamisessa tarvitaan vertailuopraattoreita.
	 x==y	yhtä suuri
	 x!=y	eri suuri
	 x>y	suurempi kuin
	 x<y	pienempi kuin
	 x>=y 	suurempi tai yhtä suuri kuin
	 x<=y 	pienempi tai yhtä suuri kuin
	 x===y 	identtisesti yhtä suuri

1 < 2 < 3

1 < 3 < 2

Loogiset operaattorit

Loogisia operaatioita tarvitaan yhdistettäessä erilaisia ehtolausekkeita. 
	Not[x] tai ! 	ei
	 x&& y  	ja
	 x || y		tai
	 Xor[x,y,...]	 joko ... tai

Not[True]

Raketin lentorata

Tarkastellaan raketin lentorataa sellaisessa tapauksessa, että raketin polttama polttoaine on merkittvä osa kokonaismassasta, joten lennon aikana raketin massan muutos täytyy huomioida.
 Massa ajan funktiona m[t_]:=If[c t<p,m0-ct,m0-p].
Toinen tapa määritellä massa ajan funktiona:

Clear[m]
m[t_] := m0 - c t /; c t < p
m[t_] := m0 - p /; c t >= p

Nousun differentiaaliyhtälö

dy := m[t] x''[t] == fm[t] - m[t] g

Työntövoima ajan funktiona

fm[t_] := f /; c t < p
fm[t_] := 0 /; c t >= p

Vakiot yksiköissä kilogramma, kilometri ja sekunti g=9.81 m/(s^2), c = 1 kg/s, m0=1000 kg, p = 100kg, f=15000N

g = 9.81*10^-3;
c = 1;
m0 = 1000;
p = 100;
f = 15;

Differentiaaliyhtälön ratkaisuna saadaan paikka=paikka(t)

rat1 = 
        NDSolve[{dy, x[0] == 0, x'[0] == 0}, 
          x[t], {t, 0, 1000}, 
          AccuracyGoal -> 15];
paikka = x[t] /. rat1[[1]];
nopeus = D[paikka, t];
Plot[paikka, {t, 0, 300}]
Plot[nopeus, {t, 80, 200}]

Lasketaan hetki ja nopeus, kun korkeus 40km

tloppu = 
      FindRoot[paikka == 40, {t, 180, 200}]
tl = t /. tloppu
loppunop = nopeus /. t -> tl

ja radan maksimi

tmax = 
      FindRoot[nopeus == 0, {t, 150, 200}]
tm = t /. tmax
xmax = paikka /. t -> tm

Polttoaineen loputtua raketti avaa laskuvarjon 40 km korkeudessa.  Laskeutumista kuvataan  differrentiaaliyhtälöllä, jossa huomioidaan myös ilmanvastus, ilmanvastuskertomeksi valitaan γ=3000 kg/s

difflasku := (m0 - p) y''[t] + 
          kitka y'[
              t] == -(m0 - p) g
kitka = 3000

Laskeutumisen ratkaisu

rat2 = 
        NDSolve[{difflasku, y[0] == 40, y'[0] == loppunop}, 
          y[t], {t, 0, 20000}];
paikka2 = y[t] /. rat2[[1]];
Plot[paikka2, {t, 0, 20}, PlotRange -> All]

Nopeus maahan tultaessa

nopeus2 = D[paikka2, t];
Plot[nopeus2, {t, 0, 10}];
tmaa = 
        FindRoot[
          paikka2 == 0, {t, 10000, 
            15000}];
tm2 = t /. tmaa
nopeusmaa = nopeus2 /. t -> tm2
kmt = 3600 nopeusmaa km/t

Piirretään koko lentorata

lentorata[tt_] := 
      If[tt ≤ tl, paikka /. t -> tt, 
        paikka2 /. t -> tt - tl]
Plot[lentorata[tt], {tt, 0, 15000}]
Plot[lentorata[tt], {tt, tl - 50, tl + 10}, 
      PlotRange -> All]

Oikea maan vetovoima

Todellisessa tilanteessa maan vetovoima ei ole vakio vaan kääntäen verrannollinen etäisyyden neliöön

<< Miscellaneous`PhysicalConstants`
rmaa = EarthRadius/Meter*10^-3;
Print[rmaa, " km"]

Gravitaatiovoima

grav := g (rmaa/(z[t] + rmaa))^2

Raketin diff.yht.

diff3 := m[t] z''[t] == fm[t] - m[t] grav

Ratkaisun ja rajanopeuden etsiminen  oikealle gravitaatiovoimalle

Jos raketin kineettinen energia, 1/2m(Subscript[v, p])^2, silloin, kun polttoaine loppuu on suurempi kuin sen potentiaalienergia, mgR^2/(R+z),  niin se vapautuu maan vetovoimasta. Rajanopeudeksi saadaan (Subscript[v, raja])=(Sqrt[2  g]) R/(Sqrt[(R + z)]).  
Polttoaineen loppumisen jälkeen raketin liikeitä kuvaa yhtälö 
z"[t]=grav, 
jonka ratkaisuna saadaan raketin nopeus (v^2)=2*g* (R^2)/(R+z) +((Subscript[v, ∞])^2). Hyvin kaukana maasta raketti matkaa siis nopeudella  (Subscript[v, ∞]). Jotta rajanopeus saavutettaisiin, täytyy työntövoimaa kasvattaa edellisestä huomattavasti ja lisätä polttoainetta rakettiin.

f = 15;
p = 650;
tpo  =  p/c;
max  = 1000;
rat3 = 
        NDSolve[{diff3, z[0] == 0, z'[0] == 0}, 
          z[t], {t, 0, max}, AccuracyGoal -> 17];
paik[t_] = z[t] /. rat3[[1]];
vraja  =  rmaa*Sqrt[2*g/(rmaa + paik[tpo])];
Print["rajanopeus = ", vraja, " km/s korkeudella " ,  
      paik[tpo], " km"]
nop[t_] = D[paik[t], t];
Print["v_p = ", nop[tpo], " km/s"];
vaareton = Sqrt[nop[tpo]^2 - vraja^2];
Print["nopeus kaukana on ", vaareton, " km/s"]
Plot[paik[t], {t, 0, max}]
Plot[nop[t], {t, 0, max}]

Tehtäviä

67. Tutustu rakettiesimerkkiin ja kokeile mitä tapahtuu, jos muutat vakioiden, kuten alussa olevan polttoaineen arvoja. 
68. Määrittel funktio f[x], joka saa arvon 1, kun x<=0 ja arvon 0, kun x>0. Piirrä funktion kuvaaja.

Switch-lause

Switch[expr, form1, value1, form2, value2, ...,_,other value] 

Switch-lause laskee lausekkeen  expr arvon ja vertaa tulosta jokaiseen formi arvoon. Jos expr=formi niin Switch-komennon tuloksena saadaan valuei. Kun mikään arvo ei käy, niin jätetään Switch-lause suorittamatta, tai jos alaviiva vaihtoehto on lisätty, niin suoritetaan alaviivan jälkeen tulevat käskyt.

Toistorakenteet

Mathematica tuntee For, Do ja While -komennot.

For-silmukka

For-komento on muotoa 
 For[alkuasetukset puolipisteillä erotettuna, ehto, laskuri, komennot puolipisteillä erotettuna].  
 Alkuasetukset suoritetaan vain kerran, ennen muita toimia. Ehto tarkistetaan joka kierroksella, ja jos se on epätosi, poistutaan silmukasta, muuten jatketaan silmukan suorittamista järjestyksessä:
	ehto
	komennot 
	laskuri,
	ehto
	jne.
Laskuri ja komento ovat tavallisia Mathematica-komentoja. Jos halutaan useampi komento, voidaan ne erottaa puolipisteellä tai kirjoittaa listaksi.

Mathematica tuntee C-kielen tyyliset lyhennysmerkinnät:
	Lyhennysmerkinnät
Lyhennys  	Merkitys
a+=b      	a=a+b
 a-=b		a=a-b
 a/=b		a=a/b
a*=b		a=a*b
++a, a++	 a=a+1
 --a,  a--	 a=a-1

Esimerkki:

For[a = 0; i = 0, i < N[Pi], i += N[Pi/5], a += Sin[i]; 
      Print[a]]

Do -silmukka

Do-komento on muotoa  
Do[komennot puolipisteillä erotettuna, laskuri, alkuarvo, loppuarvo, askel]. 
Jos  askel puuttuu käytetään oletusarvoa askel=1. Jos sekä askel että  alkuarvo puuttuvat niin käytetään molemmille oletusarvoa 1.

Esimerkki:

a = 0
Do[a += Sin[i]; Print[a], {i, 0, N[Pi], N[Pi/5]}]

While-silmukka

While-komento on muotoa  

While[ehto, komennot puolipisteillä erotettuna]

Komentojen joukossa täytyy olla lauseke, joka muuttaa ehdon totuusarvoa silmukkaa suoritettaessa.

Esimerkki:

a = 0
i = 0
While[i < N[Pi], i += N[Pi/5]; a += Sin[i]; Print[a]]

Nest ja FixedPoint

Lisäksi Mathematica osaa Nest ja FixedPoint-toistorakenteet, joihin ei tässä sen tarkemmin puututa.

Tehtäviä

69. Kuten melkein sanasta voi päätellä, eivät edellä esitetyt esimerkit tee aivan täsmälleen samaa asiaa. Mitä ne tekevät eri tavalla?
70. Tutki, mitä eroa on lausekkeilla a++ ja ++a?

Vuorovaikutus ohjelman ajon aikana

Usein on tarve pyytää käyttäjää syöttämään esimerkiksi jokin luku. Tämän voi tehdä Input-komennolla. 

Esimerkiksi

2 + Input["Anna luku! "]

Kuten aikaisemmin oli jo esillä, Mathematicassa voidaan datoja lukea myös tiedostosta ReadList-komennolla. 
Komennolla
	datat=ReadList["file", {Number,Number}]
luetaan tiedostosta  file kahteen sarakkeeseen sijoitettuja lukuja taulukkoon datat.

Tiedostoon kirjoittamista varten täytyy avata kirjoituskanava käskyllä

kanava  =  
      OpenWrite["tiedosto", FormatType -> OutputForm]

Tiedostoon kirjoitetaan kanavan kautta.

x = 2
y = 4
Write[kanava, x, "  ", y]

Lopuksi kirjoituskanava suljetaan komennolla Close:

Close[kanava]

Tehtäviä

71. Kirjoita Mathematica ohjelma, joka laskee n ensimmäistä Fibonaccin lukua. Fibonaccin luvut saadaan rekursiokaavalla
(Subscript[x, i])=(Subscript[x, i - 1])+(Subscript[x, i - 2])
(Subscript[x, 0])=0
(Subscript[x, 1])=1,
eli jokainen luku on kahden edeltäjänsä summa. Lukumäärä n kysytään käyttäjältä. Lukujen tulostukseen voit käyttää Print-komentoa.

Kertaustehtäviä

11. Määrittele paloittain funktio, joka saa arvon x^2, kun x<3, arvon x+5 kun 3<=x<7 ja arvon 14, kun x>=7. Piirrä f:n kuvaaja. Mieti, miten saisit funktiosta jatkuvan.
12. Laske summa \(TraditionalForm\`\(\(\ \)\(∑\)\)\+\(i\  = \ 0\)\%10\) i^2   käyttäen Do-ja While-silmukoita.
13. Laske rekursiokaavan
(Subscript[y, i])=3  (Subscript[y, i - 1])+4 (Subscript[y, i - 2])
(Subscript[y, 0])=1
(Subscript[y, 1])=9
mukainen (Subscript[y, i]):n arvo, kun i=7.

Schrödingerin yhtälö

Atomien ja molekyylien ominaisuuksia voidaan kuvata lineaarisella toisen kertaluvun differentiaaliyhtälöllä, jossa on ominaisarvo. Yhtälöä kutsutaan  Schrödingerin yhtälöksi. Se voidaan kirjoittaa seuraavaan muotoon 
	y''(r)- V(r) y(r)= e y(r).
Funktio V(r) kuvaa kahden mikromaailman hiukkasen ( atomi, elektroni, protoni jne) välistä vuorovaikutusta ja e  on yhtälön ominaisarvo, jonka arvo määräytyy ratkaisulle  y(r) vaadituista rajaehdoista, kun r ->∞ ja/tai r -> 0. 
Esimerkiksi vetyatomin tapauksessa Schrödingerin yhtälön ratkaisuna  saadaan aaltofunktio y(r), jonka neliö kuvaa elektronin todennäköisyysjakautumaa protonin ympäristössä, sekä ominaisarvo e  eli energia, jolla protoni sitoo elektronin itseensä.
Yhtälön  ratkaiseminen tapahtuu siten, että ensin arvataan ominaisarvo ja sitä käyttäen ratkaistaan differentialiyhtälö (NDSolve-komento). Jos ratkaisu ei toteuta haluttuja rajaehtoja, niin suoritetaan parempi arvaus ominaisarvolle ja yritetään uudelleen.

Mathematica ohjelma vetyatomin perustilalle

Edellä kuvattua yksinkertaista periaatetta voidaan kehittää siten että numeerinen tarkkuus paranee ja ratkaisusta saadaan stabiili. Tällöin pisteissä r ->0 ja r->∞ annettuja alkuarvoja käyttäen lasketut ratkaisut sovitetaan yhteen pisteessä r=(Subscript[r, v])äli siten, että sekä y(r) että sen derivaatta y'(r) ovat jatkuvia tässä pisteessä. Koska yhtälö on lineaarinen, niin ominaisarvon määräämiseen riittää, että suhde y'(r)/y(r) on jatkuva kohtauspisteessä r=(Subscript[r, v])äli. Lineaarisen differentiaaliyhtälön ratkaisunhan voi vapaasti kertoa vakiolla ja tuloksena saadaan yhtälölle uusi ratkaisu. Lopuksi ratkaisu normitetaan siten, että ∫  (y^2)(r) dr =1. 
Ominaisarvon arvausta parannetaan binaarihaun avulla. Haku toimii silloin, kun on löydetty kaksi ominaisarvoa  e=e1 ja  e=e2, jotka antavat eri merkin erotukselle 
Δ=  (Subscript[y, 0])'((Subscript[r, v])äli)/(Subscript[y, 0])((Subscript[r, v])äli)-(Subscript[y, ∞])'((Subscript[r, v])äli)/(Subscript[y, ∞])((Subscript[r, v])äli).
Merkinnät (Subscript[y, 0])((Subscript[r, v])äli)  ja  (Subscript[y, ∞])((Subscript[r, v])äli)  tarkoittavat origosta ja äärettömyydestä lähteviä ratkaisuja.

Schrödingerin yhtälö

Clear[e, y]
yhtalo :=  y''[r] - (Vcoul[r] - e) y[r] == 0

Coulombin potentiaali

Vcoul[r_] :=  -2/r  +  
        L (L + 1)/r^2

Alkuarvoissa asetetaan 
- laskentatarkkuutta kuvaava parametri, tark
- piste, josta diff. yhtälön ratkaiseminen alkaa lähellä origoa, ralku
- piste, johon diff. yhtälön ratkaiseminen päättyy, rvali
- piste, josta diff. yhtälön ratkaiseminen alkaa kohti origoa, rloppu
- ensimmäinen arvaus ominaisarvolle, e=e1, ja sitä vastaava Δ(e1)=ero1.
- toinen arvaus ominaisarvolle, e=e2, ja sitä vastaava Δ(e2)=ero2.

L = 0;
tark     =  10^(-12) ;     
ralku    =  10^(-8) ;     
rvali    =  1;            
rloppu   =  8;
e1           =  -1.1;
e2           =  -0.9;

Ratkaisufunktion asymptoottinen käyttäytyminen alkupisteessä y(r->0)=c1 * r  ja loppupisteessä y(r->∞) = c2 * (e^(-(Sqrt[-e]) r)).  Lasketaan myös funktion derivaatat alku- ja loppupisteissä.

Clear[r, e, yal, dyal, ylo, dylo]
c1  =  1
c2 = 10^-4
yal[r_]  =  c1*r
dyal[r_] = D[yal[r], r];
ylo[r_]  =  c2*E^(-Sqrt[-e]*r)
dylo[r_] = D[ylo[r], r];

Lasketaan logaritmisten derivaattojen erotus, ero1, sovituspisteessä, rvali, sen jälkeen arvataan uusi ominaisarvo, e2,  ja lasketaan sitä vastaava logaritmisten derivaattojen erotus, ero2. Mikäli ero1 ja ero2 ovat erimerkkiset, niin oikean ominaisarvon täytyy olla välillä e1 < e < e2. Parannetaan arvausta puolittamalla väli jokaisella kierroksella ja huolehtimalla siitä, että ero1 ja ero2 ovat aina erimerkkiset.

Alkuarvojen etsintä ominaisarvolle (e1 ja e2)

Clear[y, r]
e = e2
\(rata\  = \ 
      NDSolve[{yhtalo, 
          y[ralku] == yal[ralku], \(y'\)[ralku] == 
            dyal[ralku]}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ y[r], {r, ralku, 
          rvali}, AccuracyGoal -> 15]\)

kuvaalku = 
      Plot[y[r] /. rata, {r, ralku, rvali}, PlotRange -> All]

Clear[alk]
alk[r_]         =  
      y[r] /. rata[[1]]
logalk[r_]       =  D[alk[r], r]/alk[r]

Clear[lop]
ratl  =  
      NDSolve[{yhtalo, 
          y[rloppu] == ylo[rloppu], y'[rloppu] == 
            dylo[rloppu]}, y[r], {r, rloppu, rvali}, 
        AccuracyGoal -> 15]

kuvaloppu = 
      Plot[y[r] /. ratl, {r, rvali, rloppu}, PlotRange -> All]

lop[r_] =  y[r] /. ratl[[1]]
loglop[r_] =  D[lop[r], r]/lop[r]

\(\(\(\ \)\(ero\  = \ \(\(logalk[rvali] - 
        loglop[rvali]\)\(
\)\)\)\)\)

Show[kuvaalku, kuvaloppu]

Logaritmisten derivaattojen erotus arvoilla e1 (ero1) ja e2 (ero2)

ero1       =  0.174;
ero2      =  -0.197;

Varsinainen iterointi

\(
While[
      Abs[e1 - e2] > 
        tark, 
\ \ \ \ \ \ \ \ e\  = \ \((e1 + e2)\)/2; 
		 (*\ 
        alkuosa\ *) 
\ \ \ \ \ \ \ \ rata\  = \ 
        NDSolve[{yhtalo, 
            y[ralku] == yal[ralku], \(y'\)[ralku] == 
              dyal[ralku]}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ y[r], {r, ralku, 
            rvali}, AccuracyGoal -> 15]; 
\ \ \ \ \ \ \ \ alku[
          r_]\ \ \ \ \ \ \ \  = \ 
        y[r] /. rata[\([1]\)]; 
\ \ \ \ \ \ \ \ logalk[r_]\ \ \ \ \ \  = \ 
        D[alku[r], r]/alku[r]; 
		 (*\ 
        loppuosa\ *) 
\ \ \ \ \ \ \ \ ratl\  = \ 
        NDSolve[{yhtalo, 	
            y[rloppu] == ylo[rloppu], \(y'\)[rloppu] == 
              dylo[rloppu]}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ y[r], {r, 
            rloppu, rvali}, AccuracyGoal -> 15]; 
\ \ \ \ \ \ \ \ loppu[
          r_]\ \ \ \ \ \ \  = \ 
        y[r] /. ratl[\([1]\)]; 
\ \ \ \ \ \ \ \ loglop[r_]\ \ \ \ \ \  = \ 
        D[loppu[r], r]/loppu[r]; 
	\ \ \ \ \ \ \ \  (*\ 
        alku - \ 
          ja\ loppuosan\ sopivuus\ yhteen\ *) 
\ \ \ \ \ \ \ \ ero\ \ \ \ \ \ \ \ \ \ \ \ \  = \ logalk[rvali] - loglop[rvali]; 
\ \ \ \ \ \ \ \ If[
        ero\ ero1\  >= \ 0, 
\ \ \ \ \ \ \ \ e1\ \ \ \ \ \  = \ 
          e; 
\ \ \ \ \ \ \ \ ero1\ \ \ \  = \ 
          ero, 
\ \ \ \ \ \ \ \ e2\ \ \ \ \ \  = \ 
          e; 
\ \ \ \ \ \ \ \ ero2\ \ \ \  = \ 
          ero;
\ \ \ \ \ \ \ \ ];
] (*\ 
      While - silmukka\ loppuu\ *) \)
Print["ominaisarvo = ", e]    
Print["logaritmisten derivaattojen erotus = ", ero]

Sovitetaan ratkaisut yhteen ja normitetaan. Nämä ehdot määräävät vakiot c1 ja c2.

\(\ \ \ \ \ \ \ \  (*\ \ \ \ \ \ ratkaisujen\ sovittaminen\ yhteen\ \ \ \ \ \ \ *) 
\(c2\  = \ 
        alku[rvali]/
          loppu[rvali];\)
\ \ \ \ \ \ \ \  (*\ \ \ \ \ \ ratkaisun\ normittaminen\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ *) \)
\(norm\  = \ 
      NIntegrate[\((alku[r])\)^2, {r, ralku, rvali}] + 

        c2^2\ NIntegrate[\((loppu[r])\)^2, {r, rvali, rloppu}]\)

\(\(\(\ \ \ \ \ \ \ \ \)\( (*\ \ \ \ \ \ \(kuvaajat\)\(:\)\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ *) \)\(
\)\(\ \ \ \ \ \ \ \ \)\( (*\ \ \ \ \ \ Todennäköisyystiheydet\ \ \ \ \ \ \ \ \ \ *) \)\(
\)\(kuvaalku = 
      Plot[\((alku[r])\)^2/norm, {r, ralku, rvali}, 
        PlotRange -> All]\)\)\)

kuvaloppu = 
      Plot[(c2 loppu[r])^2/ norm, {r, rvali, rloppu}, 
        PlotRange -> All]

Show[kuvaalku, kuvaloppu, PlotRange -> All]

Lasketaan varsinaiset aaltofunktiot,  y(r)/r.

kuva3 = 
        Plot[(alku[r]/r)/Sqrt[norm], {r, ralku, rvali}, 
          PlotRange -> All];
kuva4 = 
        Plot[(c2 loppu[r]/r)/Sqrt[norm], {r, rvali, rloppu}, 
          PlotRange -> All];
Show[kuva3, kuva4, PlotRange -> All]

Tarkka ratkaisu

Lasketaan myös Schrödingerin yhtälön tarkka ratkaisu Coulombin potentiaalille (-2/r), kun tiedetään, että ominaisarvot ovat muotoa e=1/(k+L+1)^2, missä ns pääkvanttiluku k saa arvot  0, 1, 2 ....  Tuloksena saadaan hypergeometrisia funktioita, jotka ovat Laguerren polynomeja.

Clear[y, e, n, tarkka, L, k]
e = -1/(1 + k + L)^2;
rat =  
        DSolve[yhtalo, y[r], r]  // Simplify;
tarkka[k_, L_] = PowerExpand[rat] // Simplify

Valitaan kvanttiluvut k ja L, normitetaan ratkaisu ja piirretään kuvaaja.

k = 0;
L = 0;
fun = 
        y[r] /. (tarkka[k, L] /. 
              C[2] -> 0);
fu = fun[[1]]

normtark = Integrate[fu^2, {r, 0, Infinity}]

aalto[
        r_] = fu/(Sqrt[normtark]*
              r)//PowerExpand

Plot[aalto[r], {r, 0, 8}, PlotRange -> All]

Tehtäviä

72. Yritä löytää vetyatomin 1. viritystilan energia ja määrittää sitä vastaava aaltofunktio.

Vektoreita  ja matriiseja

Peruskomentoja

Mathematicassa vektorit ovat listoja. Esimerkiksi vektori a=(x,y,z) esitetään listana a={x,y,z}}. Matriisit vuorostaan ovat listoja, joiden alkioina on listoja. 
       			m = (({{a, b}, {c, d}}))
 annetaan Mathematicalle muodossa  m={{a,b},{c,d}}.
Tärkeimmät vektoreiden ja matriisien laskutoimitukset:
Vektorin p kertominen skalaarilla a: 	p  a 
Skalaarivakion  a lisääminen kaikkiin vektorin p komponentteihin: 	a + p
Piste- eli sisätulo:	 	a . b
Matriisin m ja vektorin a tulo:   	m . a
Kahden matriisin tulo:		m1 . m2	
Transponoitu matriisi:		Transpose[m]
Käänteismatriisi:		Inverse[m]
Determinantti:			Det[m]
Ominaisarvot:			Eigenvalues[m]
Ominaisvektorit:		Eigenvectors[m]
Matriisiyhtälön m x = a ratkaiseminen x:n suhteen:   LinearSolve[m,a]

Mathematica osaa näyttää matriisit selkeämmässä muodossa MatrixForm-komennon avulla. Tämän komennon voi kirjoittaa myös varsinaisen laskukomennon perään. Esimerkin matriisi m on määritelty kuten tehtävissä.

m = {{1, 2}, {3, 4}}

Inverse[m]  //  MatrixForm

Tehtäviä

74. Määrittele matriisi
	m=(({{1, 2}, {3, 4}}))
ja laske sen käänteismatriisi, ominaisarvot, ominaisvektorit, transponoitu matriisi ja determinantti.
75. Esitä edellisen tehtävän tulokset matriisimuodossa.
76. Määrittele vektori a=(4,5)  ja ratkaise x yhtälöstä mx=a.
77. Kokeile auttaako  MatrixForm esittämään vektoreita selkeämmässä muodossa.
78. Määrittele jokin 3x3-matriisi ja tutki, onko sillä käänteismatriisia. Jos on,niin laske se ja totea tuloksen oikeellisuus kertomalla matriisi käänteismatriisillaan.

Schrödingerin yhtälön ratkaiseminen matriisimuodosta

Esimerkkinä ratkaistaan Schrodingerin yhtälö 
	y''(r) -V(r) y(r) = e  y(r)
matriisimuodossa. Sitä varten kirjoitetaan toinen derivaatta ns. kolmen pisteen kaavan avulla,
	y''(r) -> (1/(h^2))( y(r-h) - 2 y(r) + y(r+h) )
Piste r on matriisin diagonalielementti ja pisteet r-h ja r+h vasemman ja oikeanpuoleisia sivudiagonaaleja.

Alkuarvoiksi tarvitaan
- pisteiden vali, h
- alkupiste differentiaaliyhtälön ratkaisemiselle, alku
- loppupiste differentiaaliyhtälön ratkaisemiselle, loppu
- pisteiden lukumäärä, nn

h   =  0.05;  
alku  =  0;        
loppu =  8; 
nn := (loppu - alku)/h + 1;

Potentialiksi valitaan tällä kertaa harmoonisen värähtelijän potentiaali, V(r) = (r^2).

Clear[x, r, v, i, j, d]
v := Table[r^2, {r, alku, loppu, h}];

Muodostetaan kerroinmatriisi toisen derivaatan ja potentiaalitermin  laskemiseksi.

d := -1/h^2
\(\(mat = 
        Table[Switch[i - j, \(-1\), d, 0, v[\([i]\)] - \ 2\ d, 1, d, _, 
            0]
\ \ \ \ \ \ \ \ , {i, nn}, {j, nn}];\)\)

Lasketaan ominaisarvot ja ominaisvektorit. Ratkaisuna saadut ominaisarvot on järjestetty siten, että alin ominaisarvo on viimeinen ja ylin ensimmäinen. Vastaavat ominaisvektorit ovat samassa järjestyksessä.

eigen = Eigenvalues[mat];
eigvec = Eigenvectors[mat];

Lasketaan peräkkäisten ominaisarvojen erotus.  Jos tulokset ovat tarkkoja, niin erotuksien pitäisi olla yhtäsuuria.

erotus[i_] := eigen[[i]] - eigen[[i - 1]]
ListPlot[Table[{i, erotus[i]}, {i, nn, nn - 5, -1}], 
      PlotJoined -> False]

Tulostetaan muutama viimeinen ominaisvektori siten, että interpoloidaan polynomi laskettujen pisteiden kautta.

normi = 
      NIntegrate[
        Evaluate[
          Table[Interpolation[-eigvec[[nn - i]]][(x - alku)/
                    h + 1], {i, 0, 0}]^2], {x, alku, loppu}]

kuva1 = 
      Plot[Evaluate[
          Table[Interpolation[-eigvec[[nn - i]]][(x - alku)/
                    h + 1 ]/(Sqrt[normi]), {i, 0, 0}]], {x, alku, 4}, 
        PlotRange -> All]

Tarkka ratkaisu

Schrödingerin yhtälön tarkka ratkaisu harmooniselle oskillaattorille (r^2).

Clear[n, L, harmosk, k]
e = 2 (2*k + L) + 3
harmosk[k_, 
        L_] = u''[r] + (e - r^2 - L*(L + 1)/r^2) u[r] == 
        0
rathar[k_, L_] = DSolve[harmosk[k, L], u[r], r]

Clear[aahar]
k = 0
L = 0
rat = 
      rathar[k, L] /. C[2] -> 0 // Simplify
aah = u[r] /. rat[[1]]
no = Integrate[aah^2, {r, 0, Infinity}]
aalhar[r_] = 
      aah/(Sqrt[no]) // PowerExpand
kuva2 = Plot[aalhar[r], {r, 0, 4}, PlotRange -> All]

Show[kuva1, kuva2]