2008-11-19

5. Fizika i Mathematica

5.1 Klasična mehanika

Gibanje u centralnom sferno-simetricnom potencijalu

Zadatak 5-1:    Slično kao na kraju drugog poglavlja, napravite simulaciju gibanja čestice u centralnom sferno-simetričnom gravitacijskom potencijalu

V(r)=-(1/r)

Igrajte se s početnim uvjetima da dobijete: (a) zatvorenu (eliptičnu) putanju (b) otvorenu (hiperboličnu) putanju.

Naputak: Kako znamo da je gibanje u ravnini,  dovoljno je raditi 2D simulaciju tj. uzeti r=(Sqrt[(x^2)+(y^2)]).

Solitoni

Solitoni su stabilna lokalizirana rjesenja nelinearnih diferencijalnih jednadzbi i javljaju se u mnogim podrucjima fizike. Kao  primjer uzmimo tzv. Korteweg de Vriesovu (KdV) jednadzbu

((∂u)/(∂t))-6 u ((∂u)/(∂x))+(((∂^3)u)/(∂(x^3)))=0

Mathematica moze numericki rjesavati neke parcijalne nelinearne jednadzbe, ali za pronalazenje solitonskih rjesenja potrebne su specijalne tehnike u koje necemo ovdje ulaziti, vec cemo samo provjeriti i graficki prikazati poznata rjesenja iz literature

Zgodno je definirati diferencijalni operator koji djeluje na funkciju u(x,t):

KdV[u_]:=D[u,t]-6 u D[u,x]+D[u,x,x,x]

Pokusajmo s rjesenjem oblika:

utry[x_,t_]:=a ((Cosh[x-b t])^(-2))

KdV[utry[x,t]] //Simplify

Da bi to bilo nula (dakle rjesenje KdV jedandzbe) za svaki x i t, ocito mora biti b=4 ...

%/.b->4

... i a=-2.

%/.a->-2

Dakle rjesenje je:

u1[x_,t_]:=-2 ((Cosh[x-4 t])^(-2))

KdV[u1[x,t]]//Simplify

Zadatak 5-2: Zgodno je prikazati oblik rjesenja u razna vremena, i animirati dobivene grafove. Upotrijebite 
funkciju Animate, da graficki prikazete soliton u1[x,t], za x izmedju -20 i 20, te za vrijeme t izmedju -4 i 4.

Zadatak 5-3: Uvjerite se da je dolje definirana funkcija u2[x,t] takodjer soliton koji zadovoljava KdV jednadzbu. Animirajte je isto kao gore.

Zadatak 5-4: Tzv. sinus-Gordonova diferencijalna jednadzba je

(((∂^2)u)/(∂(t^2)))-(((∂^2)u)/(∂(x^2)))+sin(u)=0

Uvjerite se da je tzv. "kink"

4 arctan( exp(((x-v t)/(Sqrt[1-(v^2)]))))

rjesenje ove jednadzbe (Naputak: FullSimplify). Animirajte ga.

5.2 Elektrodinamika

Potencijal staticke raspodjele naboja

Zadatak 5-5:  Potencijal ϕ  u tocki (Overscript[r, →])=(x, y, z) koji je posljedica staticke raspodjele N nabijenih cestica s nabojima (Subscript[q, 1]), (Subscript[q, 2]),..., (Subscript[q, N]) i polozajima (Subscript[(Overscript[r, →]), 1])=((Subscript[x, 1]), (Subscript[y, 1]), (Subscript[z, 1])),  (Subscript[(Overscript[r, →]), 2])=((Subscript[x, 2]), (Subscript[y, 2]), (Subscript[z, 2])), ..., (Subscript[(Overscript[r, →]), N])=((Subscript[x, N]), (Subscript[y, N]), (Subscript[z, N]))   je dan formulom

ϕ(x,y,z) = (∑_(n=1)^N)((Subscript[q, n])/(|(Overscript[r, →])-(Subscript[(Overscript[r, →]), n])|))

Definirajte funkciju plotPotencijal[naboji_,xgranice_,ygranice_]:=... koja crta ekvipotencijalne konture u ravnini z=0,  gdje su naboji zadani kao lista {{q1,x1,y1,z1},{q2,x2,y2,z2}, ...},

Naputak: Koristite ContourPlot da nacrtate funkciju ϕ(x,y,0).

Gibanje nabijene čestice u homogenom magnetskom polju

Zadatak 5-6: Jednadzba gibanja nabijene čestice u elektromagnetskom polju je

m (Overscript[a, -> ]) = q ((Overscript[E, ->])+(Overscript[v, ->])× (Overscript[B, ->]))

Riješite je za slučaj homogenog magnetskog polja duz z-osi i nacrtajte putanju za neke tipične početne uvjete.

Naputak: Mozete postupiti ovako:

1. Definirajte vektore brzine, akceleracije i magnetskog polja

brzina={x'[t],y'[t],z'[t]};
akceleracija={x''[t],y''[t],z''[t]};
magnetskoPolje={0,0,B};

Tada su diferencijalne jednadzbe gibanja

jednadzbeGibanja= m akceleracija == q Cross[brzina,magnetskoPolje]//Thread

(Pogledajte u manualu upotrebu funkcije  Thread.Koja je funkcija gore  "threaded"?)

2. Početne uvjete definirajte kao posebnu listu jednadzbi.

3. Riješite diferencijalne jednadzbe upotrebom DSolve

4. Nacrtajte rješenje upotrebom ParametricPlot3D

Zadatak 5-7:    Upotrebom paketa Graphics`PlotField3D` superponirajte na gornju putanju skicu vektorskog magnetskog polja

5.3 Kvantna fizika

Bohrov model atoma

Zadatak 5-8: Koristeći Bohrov uvjet za orbitu elektrona u atomu vodika

moment impulsa = n 

izračunajte i nacrtajte spektar vodika.

Naputak: Možete postupiti ovako:

Gornji uvjet

eq1=m v r == n 

i jednadžba gibanja (Coulombova sila = m*centripetalno ubrzanje)

eq2=((e^2)/(4 π (Subscript[ϵ, 0])(r^2)))==m ((v^2)/r)

daju dvije jednadžbe s dvije nepoznanice koje se mogu riješiti po r i v. Ta riješenja se ona mogu uvrstiti u izraz za ukupnu energiju E(n)=kinetička + potencijalna, što daje Bohrov spektar. Pomoću paketa PhysicalConstants`  i Units` treba napraviti konverziju konstanti da se dobije energija u elektron-voltima (eV).

5.4 Astronomija

Efemeride

Zadatak 5-10: U datoteci bsc.dat, koja se moze skunuti s http://www.phy.hr/~kkumer/mma/files.tgz  nalaze se koordinate oko 9000 najsjajnijih zvijezda u formatu:

[redni broj] [sati RA] [minute RA] [sekunde RA] [stupnjevi dec] [minute dec] [sekunde dec] [magnituda]

gdje je RA= rektascenzija ("right ascension"), a dec=deklinacija (declination). Magnituda opisuje sjaj zvijezde (veca magnituda=manji sjaj). Npr. redak pod rednim brojem 2491:

2491   06 45 08.9     -16 42 58     -1.46

kaze da se na koordinati s rektascenzijom 6 sati 45 min 8.9 sec i deklinacijom -16 stupnjeva 42 minute i 58 sekundi nalazi zvijezda s magnitudom -1.46.  To je, usput budi receno, najsjanija zvijezda u katalogu, Sirius.

Ucitajte katalog u neku listu i nacrtajte kartu cijelog neba gdje ce sjajnije zvijezde biti predstavljene vecim tockicama.  Omogucite zadavanje granicne magnitude (tako da se iscrtavaju samo zvijezde sjajnije od neke granice) i crtanje samo dijela neba. Prepoznajte neko sazvjezdje. Nemojte pokusavati crtati 3D "globus", vec 2D kartu!

Za nadogradnju, mozete pokusati (1) napraviti funkciju koja daje izgled neba s neke geografske sirine za neko doba godine ili pak (2) crta putanju nekog planeta.

Neki linkovi:
 http://hr.wikipedia.org/wiki/Nebeski_ekvatorski_koordinatni_sustav
 http://hpd.botanic.hr/ast/astronomija/polozaj/koor_sus.html
 http://www.physics.csbsju.edu/astro/sky/sky.ra.dec.html

Otkrivanje ekstrasolarnih planeta (spektralna analiza)

Jedna od metoda za otkrivanje planeta je opazanje promjena na zvijezdi oko koje oni orbitiraju, a koje mogu nastati bilo tako da planet zasjeni zvijezdu (i izazove promjene sjaja), bilo tako da zbog gravitacijskog utjecaja planeta zvijezda orbitira ne oko svojeg sredista vec oko zajednickog centra masa. U ovom zadatku analizirat cemo izmjerene pulseve s jednog pulsara. Pulsari obicno salju izvanredno pravilne snazne elektromagnetske pulseve, ali ukoliko oko pulsara orbitira planet, u periodu pulseva nastaju varijacije.

Nas pulsar ima osnovni period pulsiranja 6.22 milisekunde, ali nizom mjerenja u godinama 1990-1993. dobivena su odstupanja zabiljezena u datoteci dostupnoj na adresi  ftp://ftp.physics.uwa.edu.au/pub/Data/pulsar.dat

(1) Skinite datoteku i ucitajte podatke u listu pulsar. Zanimljivi su nam drugi stupac koji predstavlja vrijeme mjerenja izrazeno u decimalnim godinama (npr. 1991.3 je  (3/365)-ti dan 1991. godine), i treci stupac koji predstavlja odstupanje od osnovnog perioda pulsacije u milisekundama. Kreirajte listu varijacije={{(Subscript[t, 1]),f((Subscript[t, 1]))},{(Subscript[t, 2]),f((Subscript[t, 2]))},...} samo sa ta dva stupca i to tako da u  stupcu s vremenima mjerenja oduzmete 1990 godina radi bolje preciznosti kasnijeg racuna. Pokusajte, ako mozete, ne koristiti petlje vec funkcionalne metode ili transformacijska pravila. Prikazite ove podatke graficki.

Ocekuje se da su perturbacije uzrokovane planetima periodicke tako da varijacije mozemo opisati superpozicijom jednostavnih funkcija oblika  f(t)=a sin(ω t) + b cos(ω t), gdje je ω frekvencija orbitiranja planeta koji proizvodi perturbacije. Obzirom da planete tek trebamo "otkriti" ne znamo frekvencije. One se mogu odrediti Fourierovom analizom podataka iz liste varijacije, ali i nesto jednostavnijom metodom koju cemo sad predstaviti:

(2) Nacinite funkciju ft[ω_] koja za pretpostavljeni zadani ω, fita (funkcija Fit[])  gore spomenutu funkciju f(t) na podatke tj. odredjuje koeficijente a i b.

Takozvana spektralna energija je dana funkcijom (Subscript[E, spektralna])(ω)=(1/2)(((a(ω))^2)+((b(ω))^2)) i ona pokazuje izrazene vrhove na frekvencijama ω koje odgovaraju periodima T orbitiranja perturbirajucih planeta: T=2 π/ω.

(3) Definirajte funkciju spektar[ω_] i nacrtajte spektralnu energiju kao funkciju od ω, za ω od 1 do 100, identificirajte frekvencije {(Subscript[ω, 1]),(Subscript[ω, 2]),...} istaknutih vrhova te tako odredite broj perturbirajucih planeta. Odredite i njihove periode {(Subscript[T, 1]), (Subscript[T, 2]),...}. Poziciju vrhova mozete odrediti od oka, ali pokusajte i smisliti neku sofisticiraniju metodu.

Otkrivanje ekstrasolarnih planeta (diskretni Fourierov transformat)

Rijesit cemo identican problem, ali ovaj put metodama Fourierove analize. Poznato je da kad je neka funkcija vremena f(t) priblizno jednaka superpoziciji periodicnih funkcija (npr. sin(ω t),  exp(i ω t), ...) razlicitih perioda T=2 π/ω, onda njen  Fourierov transformat F(ω) (definiran dolje), koji je funkcija frekvencije, ima izrazene vrhove na frekvencijama ω koje odgovaraju tim periodima (vidi dolje za tocnu korespondenciju). Tako, buduci da smo ustanovili da se varijacije perioda pulsara, zbog utjecaja planeta, mogu opisati superpozicijom  periodicnih funkcija, s periodima dobivenim u prvom dijelu zadataka, ocekujemo da cemo iz Fourierovog transformata
tih varijacija opet moci ocitati frekvencije ophodnje tih planeta, slicno kao ranije iz spektralne energije.

Radit cemo tzv. metodom diskretnog Fourierovog transformata, koji je definiran kao lista {F((Subscript[ω, 1])), F((Subscript[ω, 2])), ... F((Subscript[ω, n]))} s elementima danim formulom:

F((Subscript[ω, s])) = (1/(Sqrt[n]))(∑_(r=1)^n)f((Subscript[t, r])) (e^(i (Subscript[ω, s])(r-1)));      (Subscript[ω, s])=((2π (s-1))/n)

gdje su f((Subscript[t, r])) vrijednosti funkcije f(t) uzorkovane u ekvidistantnim vremenskim trenucima {(Subscript[t, 1]), (Subscript[t, 2]),...,(Subscript[t, n])}.

(1) Isprogramirajte funkciju MyFourier[list_] koja kao argument prima listu {f((Subscript[t, 1])), f((Subscript[t, 2])), ... f((Subscript[t, n]))}, a vraca odgovarajuci diskretni Fourierov transformat {F((Subscript[ω, 1])), F((Subscript[ω, 2])), ... F((Subscript[ω, n]))}. Testirajte je prvo na jednostavnim funkcijama poznatih perioda i uvjerite se da Fourierovi transformati (nacrtani pomocu ListPlot) pokazuju izrazene vrhove na frekvencijama

ω =(( 2 π)/(ukupno vrijeme uzorkovanja))(polozaj vrha-1)

Ovdje je "polozaj vrha" jednak cijelom broju k,  tako da je F((Subscript[ω, k])) polozaj maksimuma Fourierovog transformata.
Usporedite Vasu funkciju s ugradjenom funkcijom Fourier koja radi to isto, vidi slijedeci primjer. Kod tog testiranja je potrebno prvo uzorkovati funkciju koju analizirate. Npr. uzmimo funkciju

f1[t_]:=2.1Sin[7 t]+0.8 Sin[11.5 t]+1.5 Sin[13.2 t]

Uzimamo uzorak na ekvidistantnim vremenskim razmacima:

uzorak=Table[f1[t],{t,1,8,0.1`}];
ListPlot[uzorak,Joined->True]

Transformat je onda dan kao lista (uzimamo apsolutne vrijednosti jer transformat sacinjavaju kompleksni brojevi, ali nas informacija o fazi ovdje ne zanima):

transformat=Abs[Fourier[uzorak]];
ListPlot[transformat,Joined->True,PlotRange->All]

Zbog matematickih svojstava diskretnog Fourierovog transformata, on je simetrican i informacija koja nas zanima se nalazi samo u jednoj polovici transformata, gdje vidimo tocno tri vrha, s visinama koje odgovaraju amplitudama harmonika funkcije f1[t]. Da bismo tocnije odredili pozicije vrhova trebamo povecati taj dio slike. To cemo napraviti tako da napravimo funkciju koja interpolira tocke liste i crtamo nju.

interp=Interpolation[transformat]

Plot[interp[ω],{ω,5,18},PlotRange->All]

Vidimo da su vrhovi na pozicijama 9, 14 i 16. To prema gornjoj formuli odgovara frekvencijama (primijetite da smo uzorak funkcije gore uzimali od prve do osme vremenske jedinice (nismo se odredili da li su to sekunde, godine, ...).

((2 π)/(8-1))({9,14,16}-1.)

Ovo je solidna aproksimacija tocnih frekvencija koje su (vidi gore) 7, 11.5 i 13.2  (vremenskih ((jedinica))^(-1)).

Napomena: Preciznost metode ovisi o vremenu uzorkovanja: sto duze uzimamo uzorke moci cemo bolje odrediti frekvenciju. S druge strane frekvencija uzorkovanja (razmak tocaka u uzorku) ogranicava frekvencije koje mozemo
vidjeti. Prema tzv. Nyquist-Shannon teoremu frekvencija uzorkovanja mora biti barem dva puta veca od maksimalne frekvencije signala koji uzorkujemo. (Zbog toga se kod digitalizacije muzike  koristi frekvencija uzorkovanja od 44 kHz, koja je otprilike duplo veca od maksimalne frekvencije koju cuje ljudsko uho.)

(2) Kad ste se uvjerili da Vasa funkcija MyFourier radi dobro primijenite je na analizu podataka iz liste varijacije iz prethodnog zadatka. Kako vremenski trenuci u toj listi nisu ekvidistantni, morate prvo interpolirati kontinuiranu funkciju kroz te podatke, koristeci Interpolation, a zatim uzorkovati tu funkciju. Uvjerite se da dobijete iste frekvencije i periode planeta (izrazene u danima!) kao u proslom dijelu zadatka.

(3) Kao bonus, mozete pokusati  nekako odrediti ocekivanu gresku dobivenih perioda.