DY:n ratkaiseminen osa 2

# yhtaloryhmat

Yhtälöryhmistä

Joskus ratkaistavana on useita toisistaan riippuvia funktioita, joita sitoo useampi DY, eli DY-ryhmä.

Myös yksi korkeamman asteen DY on toisinaan hyödyllisempää esittää monen ensimmäisen asteen DY:n muodostamana ryhmänä.

Esimerkki 1: DY-pari

Ratkaistaan DY-pari $\begin{cases} \green{x'}(\blue t)=\red y(\blue t)+4\green x(\blue t)\\ \red {y'}(\blue t)=\red y(\blue t)-2\green x(\blue t). \end{cases}$ Ylemmästä yhtälöstä saadaan $\red y(\blue t)=\green{x'}(\blue t)-4\green x( \blue t),$ josta voidaan ratkaista derivaatta $\red {y'}(\blue t)=\green{x''}(\blue t)-4\green{x'}(\blue t).$ Kun sijoitetaan nämä alempaan yhtälöön saadaan $\green{x''}(\blue t)-4\green{x'}(\blue t)=\green{x'}(\blue t)-4\green x(\blue t)-2\green x(\blue t),$ eli VLHDY $\green{x''}(\blue t)-5\green{x'}(\blue t)+6\green x(\blue t)=0.$ Tämän ratkaisu on (katso osa 1) $\green x(\blue t)=C_1e^{3\blue t}+C_2e^{2\blue t}$ ja sen derivaatta $\green {x'}(\blue t)=3C_1e^{3\blue t}+2C_2e^{2\blue t}.$ Sijoittamalla saadaan ratkaisu myös funktiolle $\red y(\blue t).$ DY-parin ratkaisu on yhtälöpari $\begin{cases} \green x(\blue t)=C_1e^{3\blue t}+C_2e^{2\blue t}\\ \red y(\blue t)=-C_1e^{\blue t}-2C_2e^{2\blue t}. \end{cases}$

Esimerkki 2: Korkeamman asteen DY:n muuttaminen DY-ryhmäksi

Muutetaan 2. kertaluvun DY $\red y''(\green x)=\sin(\green x)\red y'(\green x)+\green x^2\red y(\green x)$ DY-ryhmäksi.

Merkitään $\blue{a}(\green x)=\red y(\green x)$ ja $\violet b(\green x)=\red y'(\green x)$ . Näin saamme DY-ryhmän $\begin{cases} \blue{a'}(\green x)=\violet b(\green x)\\ \violet {b'}(\green x)= \sin(\green x)\violet b(\green x)+\green x^2\blue{a}(\green x), \end{cases}$ jossa esiintyy vain 1. kertaluvun DY:itä.

Tämä voidaan ratkaista numeerisesti. Ratkaisut voidaan myös esittää graafisesti $(\blue{a}(\green x),\violet b(\green x))$ -koordinaatistossa eli faasikaaviossa.

Sovellus: Reaktioyhtälö

Reaktiossa kaksi ainetta muuttuvat toisikseen nopeuksilla k_1 ja k_2, siis $\text{\red A} \overset{k_1}{\underset{k_2}{\rightleftharpoons}}\text{\blue B}.$ Funktiot $\red {\text{[A]}}(\green t)$ ja $\blue {\text{[B]}}(\green t)$ kuvastavat aineiden konsentraatioita ajanhetkellä $\green t.$ Nämä toteuttavat DY:n $\frac{\text{d}\red {\text{[A]}}}{\text{d}t}=-k_1\red {\text{[A]}}(\green t)+k_2\blue {\text{[B]}}(\green t).$ Sanallisesti: Negatiivinen termi tarkoittaa, että aineen $\text{\red A}$ häviämisnopeus kasvaa lineaarisesti konsentraation mukana. Posiitiven termi puolestaan tarkoittaa, että aineen $\text{\red A}$ muodostumisnopeus kasvaa konsentraation $\text{\blue {[B]}}$ mukana. Nämä nopeudet riippuvat myös reaktionopeuksista k_1 ja k_2.

Tässä reaktiossa aineen kokonaiskonsentraatio ei muutu, sillä massan säilymislaki on $\red {\text{[A]}}(\green t)+\blue {\text{[B]}}(\green t)=\red {\text{[A]}}(\green 0)+\blue {\text{[B]}}(\green 0)$ kaikilla $\green t$ .

Tästä saadaan DY rajoitteineen: $\begin{cases} \red {\text{[A]}}'(\green t)=-k_1\red {\text{[A]}}(\green t)+k_2\blue {\text{[B]}}(\green t)\\ \red {\text{[A]}}(\green t)+\blue {\text{[B]}}(\green t)=\violet{\text{[A]}_0}+\violet{\text{[B]}_0}\\ \red {\text{[A]}}(\green0)=\violet{\text{[A]}_0}\\ \blue {\text{[B]}}(\green 0)=\violet{\text{[B]}_0}. \end{cases}$

Sijoittamalla toisesta rivistä ratkaistu $\blue{\text{[B]}}(\green t)$ ensimmäiseen saadaan DY $\red {\text{[A]}}'(\green t)=-k_1\red {\text{[A]}}(\green t)+k_2(\violet{\text{[A]}_0}+\violet{\text{[B]}_0}-\red {\text{[A]}}(\green t)),$ eli $\red {\text{[A]}}'(\green t)=-(k_1+k_2)\red {\text{[A]}}(\green t)+k_2(\violet{\text{[A]}_0}+\violet{\text{[B]}_0}).$ Tämä on ensimmäisen kertaluvun lineaarinen epähomogeeninen DY, jonka ratkaisuksivertaa vakion variointi erimerkki 1 saadaan $\red {\text{[A]}}(\green t)=\blue Ce^{-\green t(k_1+k_2)}+\frac{k_2(\violet{\text{[A]}_0}+\violet{\text{[B]}_0})}{k_1+k_2}.$ Tästä nähdään, että konsentraation aikariippuvuus on eksponentiaalista! Tämän ei pitäisi yllättää, sillä jo alussa tiedettiin, että aineen $\text{\red A}$ häviämisnopeus riippui sen konsentraaatiosta.

Sijoittamalla tämä DY-ryhmän toiseen yhtälöön saadaan $\blue {\text{[B]}}(\green t)=\violet{\text{[A]}_0}+\violet{\text{[B]}_0}-\blue Ce^{-\green t(k_1+k_2)}-\frac{k_2(\violet{\text{[A]}_0}+\violet{\text{[B]}_0})}{k_1+k_2}.$ Myös aineen $\text{\blue B}$ konsentraatiossa on aika eksponentissa, mutta termin etumerkki on eri!

Tästä voidaan laventaa $\violet{\text{[A]}_0}+\violet{\text{[B]}_0}$ tekijällä k_1+k_2, eli $\blue {\text{[B]}}(\green t)=-\blue Ce^{-\green t(k_1+k_2)}+\frac{k_1(\violet{\text{[A]}_0}+\violet{\text{[B]}_0})}{k_1+k_2}.$ Ratkaistaan vielä $\blue C$ . Tähän riittää vain toinen yhtälö, esimerkiksi $\red {\text{[A]}}(\green t)$ : $\red {\text{[A]}}(\green 0)=\blue Ce^{-\green0(k_1+k_2)}+\frac{k_2(\violet{\text{[A]}_0}+\violet{\text{[B]}_0})}{k_1+k_2}=\violet{\text{[A]}_0}.$ $e^{\green 0}=1,$ joten $\violet{\text{[A]}_0} =\blue C+\frac{k_2(\violet{\text{[A]}_0}+\violet{\text{[B]}_0})}{k_1+k_2}$ ja vakio on $\blue C=\frac{k_1\violet{\text{[A]}_0}-k_2\violet{\text{[B]}_0}}{k_1+k_2}.$ DY-ryhmän ratkaisu on $\begin{cases} \red {\text{[A]}}(\green t)=\blue Ce^{-\green t(k_1+k_2)}+\dfrac{k_2(\violet{\text{[A]}_0}+\violet{\text{[B]}_0})}{k_1+k_2}\\ \blue {\text{[B]}}(\green t)=-\blue Ce^{-\green t(k_1+k_2)}+\dfrac{k_1(\violet{\text{[A]}_0}+\violet{\text{[B]}_0})}{k_1+k_2}\\ \blue C=\dfrac{k_1\violet{\text{[A]}_0}-k_2\violet{\text{[B]}_0}}{k_1+k_2}. \end{cases}$ Tästä voidaan osoittaa, että reaktion tasapainovakio $K_{eq}$ eli $\dfrac{\blue {\text{[B]}}(\green t)}{\red {\text{[A]}}(\green t)},\quad\text{kun}\quad \green t \to \infty$ on $K_{eq} = \dfrac{k_1}{k_2}.$

Reaktioyhtälö visuaalisesti

Tämä geogebra piirtää reaktioyhtälön ratkaisut $\red{[A]}(\green t)$ ja $\blue{[B]}(\green t).$ Punainen käyrä on siis konsentraatio $\red{[A]}$ ja sininen $\blue{[B]}.$ Vaaka-akselilla on aika. Voit muuttaa alkukonsentraatioita raahaamalla sinistä tai punaista pistettä pystyakselilla sekä säätää nopeusvakioita ja k2.

# geoReaktioYhtalo

Sovellus: Heiluri ja faasikaavio

Heilurin liikeyhtälö ilman pienen kulman aproksimaatiota on $\red{\theta''}(\green t)+\frac{g}{l}\sin(\red\theta(\green t))=0.$ Tämä voidaan esittää myös yhtälöryhmänä asettamalla $\blue v(\green t)=\red {\theta'}(\green t)$ $\begin{cases} \blue{v'}(\green t)+\dfrac{g}{l}\sin(\red \theta(\green t))=0\\ \red \theta'(\green t)=\blue v(\green t). \end{cases}$ Tätä ei voi ratkaista analyyttisesti, mutta voimme kuitenkin hahmotella yhtälön määräämiä ratoja $(\blue v,\red \theta)$ faasikaaviossa.

Yhtälöryhmästä saadaan $\frac{\blue {\blue{v'}}(\green t)}{\red {\theta'}(\green t)} =\frac{\text{d}\blue v/\text{d}\green t}{\text{d}\red \theta/\text{d}\green t} =\frac{\text{d}\blue v}{\text{d}\red \theta}=\frac{-g\sin(\red \theta)}{l\blue v}.$ Huom: "supistettaessa" pois $\text d \green t$ aikariippuvuus ja samalla nopeuden etumerkki katoaa. Enää ei siis tarkalleen tiedetä, miten systeemi käyttäytyy ajan suhteen.

Separoimalla saadaan ratayhtälö $(\blue v,\red \theta)$ $\frac{l}{2}\blue v^2=g\cos(\red \theta)+C.$ Kirjoittamalla kulmataajuus $w_0=\sqrt{\dfrac{g}{l}}$ saadaan tämä muotoon $\blue v^2=2w_0^2\cos(\red \theta)+C.$

Tämä Geogebra havainnollistaa heiluria ja sen faasikaavioita, eli ratoja faasiavaruudessa, kun w_0^2=1. Tässä faasiavaruudessa pystyakseli kuvastaa heilurin nopeutta $\blue v$ ja vaaka-akseli kiertokulmaa $\red \theta.$

"Näytä heiluri" valitsimesta voit vaihtaa heiluri- ja faasikaavionäkymän välillä. Voit zoomata ja liikuttaa koordinaatistoa hiirellä.

Heilurinäkymässä

Siirtämällä pistettä $A = (\red \theta, \blue v)$ voit tarkastella eri snapshotteja ( $\green t$ on vakio) heilurin liikeestä. Liikuttamalla pistettä rataa pitkin myötäpäivään näet, miten heiluri liikkuu ajan kuluessa. Piste siis liikkuu vaaka-akselin yläpuolella vasemmalta oikealle ja alapuolella oikealta vasemmalle. Pisteen liike on rajoitettu pisteen määrämälle radalle, koska se määrittää alkuehdon $(\red \theta_0, \blue v_0)$ heilurille.

Heilurista lähtevä nuoli kuvaa heilurin nopeutta $\blue v,$ eli samaa kuin pisteen pystyakselin koordinaatti. Huom: heiluri (varsi, pallo, nuoli) ei "elä" samassa koordinaatistossa radan kanssa, vaan on yksinkertaisesti heiluri kuvattuna sivustapäin.

Faasikaavionäkymässä

Pistettä siirtämällä voit tarkastella heilurin faasikaavioiden joukkoa, eli mitä erilaisia kaavioita syntyy, kun alkuehtoa $B=(\red \theta_0, \blue v_0)$ muutetaan.

Suljetuissa käyrissä heiluri heiluu edestakaisin: vaaka-akselin koordinaatti eli kiertokulma $\red \theta$ kasvaa ja pienenee jaksollisesti.
Muissa käyrissä heiluri kiertää täysiä kierroksia: nopeus $\blue v$ on aina positiivinen tai negatiivinen, joten kiertokulma $\red \theta$ vain kasvaa (tai pienenee) ajan kuluessa.

Minkälainen kaavio syntyisi pienen kulman approksimaatiolla $\frac{g}{l}\sin(\red\theta(\green t))\approx\frac gl\red\theta(\green t)?$

# geoVaihekuva

# numeerinen

Numeerisia ratkaisumenetelmiä

Numeerinen menetelmä alkuarvotehtäville

Eräs yksinkertaisimmista tavoista ratkaista DY numeerisesti on eulerin metodi. Kirjoittamalla DY ensimmäisen kertaluvun DY-ryhmänä voimme soveltaa eulerin metodia myös korkeamman asteen DY:ille.

Menetelmän idea

Funktion derivaatta pisteessä x_0 kertoo tähän pisteeseen piirretyn tangenttisuoran kulmakertoimen. Itse tangenttisuora funktiolle saadaan yhtälöstä y(x)=f(x_0)+f'(x_0)(x-x_0). Tämä tangettisuora antaa hyvän approksimaation funktiolle pisteen x_0 lähellä. Näin voimme arvioda funktiota numeerisesti siirtymällä pienen askeleen tämän tangenttisuoran suuntaan ja piirtämällä sinne uuden tangettisuoran. Saamme rekursiokaavan $y((n+1)\delta)\approx y(n\delta)+\delta y'(n\delta),$ jonka avulla voimme arvioida funktiota murtoviivalla, kunhan tiedämme lähtöarvon ja funktion derivaatan. Tätä voidaan edelleen soveltaa ensimmäisen kertaluvun yhtälöryhmiin.

Alla olevan animaation saat käyntiin painamalla vasemmassa alanurkassa olevaa play-nappia. Kuvaan piirtyy sinikäyrä ja sen tangenttisuora. Approksimaatio kulkee tangenttisuoran suntaan yhden askeleen dt verran ja muuttaa tämän jälkeen suuntaansa senhetkiseen sijaintiinsa piirrety tangenttivektorin mukaan. Voit myös liikuttaa pistettä A ja muuttaa askelkokoa dt.

# geoEulerSelitys

Sovellus: Toteutus pythonilla

Alla olevaa python-koodia muokkaamalla voit piirtää kahden ensimmäisen asteen DY:n muodostaman ryhmän $\begin{cases} \red y'(t)=f(\green x(t),\red y(t),t)\\ \green x'(t)=g(\green x(t),\red y(t),t) \end{cases}$ kuvaajia. Painamalla "Näytä koko koodi" voit tarkastella mitä koodissa tapahtuu, "#" merkillä alkavat rivit eivät ole varsinaista koodia vaan kommentteja, joissa selitetään mitä jokin rivi tekee.

Funktion $\green x(t)$ -kuvaaja on vihreä ja funktion $\red y(t)$ punainen. Faasikaavio piirtyy $(\green x,\red y)$ -tasoon sinisellä.

# dykuvaaja

Mahdollisia virheilmoituksia:

No file: ... Laskeminen vei liikaa aikaa, pienennä askelmäärää.

RuntimeWarning: overflow... Luvut kasvoivat liian suuriksi, pienennä askelmäärää tai alkuarvoja ja tarkista, että kirjoitit lausekkeet oikein.

Heiluri numeerisesti

Voit tarkastella heilurin liikeyhtälön ratkaisuja kopioimalla tämä yllä olevaan työvälineeseen.

# Derivaattafunktiot, kirjoita haluamasi funktio returnin perään.
# lisäämällä eteen "m." voit käyttää useimpia matemaattisia
# funktioita esimerkiksi m.sin(x)

def Dy(y,x,t):
    return -m.sin(x)

def Dx(x,y,t):
    return y


# Alkuarvot ajanhetkellä t=0
y_0 = 0
x_0 = 3


# Näytetäänkö tämä kuvaaja? True=Kyllä, False=Ei, huomaa iso alkukirjain!
faasikaavio           = False
funktion_x_kuvaaja    = True
funktion_y_kuvaaja    = True

# Askelkoko, mitä pienempi sitä tarkempi kuvaaja
dt          = 0.0005

# Askelien määrä, mitä suurempi, sitä pidempi väli kuvataan
askelia     = 60000

Funktion $\green x(t)$ kuvaaja kertoo heilurin kiertokulman ajan suhteen. Vaihtamalla $\red y'=-\green x$ saat vastaavan kuvaajan pienen kulman approksimaatiolla, joka menee tässä tilanteessa pahasti pieleen.

Vastaavan kuvaajan saat myös määrittämällä alkutilanteeksi $y(0)=\sqrt2\sqrt{1-\cos(3)}\approx1.995$ ja x(0)=0. Tämä seuraa mekaanisen energian säilymislaista.

Vaihda $\text dt\text:$ n arvoksi 0.001 miksi kuvaaja ei enää näytä järkevältä?

Aseta y(0)=2 , x(0)=0 dt= 0.00004 ja askelia= 200000. Mikä on tämän tilanteen fysikaalinen tulkinta?

Pakotettu värähtelijä numeerisesti

Tällä geogebralla voit piirtää pakotetun värähtelijän analyyttisesti saadut ratkaisut. Numeerisilla menetelmillä pääsee saman näköisiin tuloksiin pienemmällä vaivalla. Liikeyhtälö tässä tilanteessa on $x''(t)+2bx'(t)+8^2x(t)=F\cos(Tt),$ jota vastaa DY-pari $\begin{cases} y'(t)=-2by(t)-8^2x(t)+F\cos(Tt)\\ x'(t)=y(t). \end{cases}$ Voit kopioida seuraavan pätkän, joka vastaa geogebran tilannetta T=3, b=0.2 ja F=20.

# Derivaattafunktiot, kirjoita haluamasi funktio returnin perään.
# lisäämällä eteen "m." voit käyttää useimpia matemaattisia
# funktioita esimerkiksi m.sin(x)

def Dy(y,x,t):
    return -64*x-0.4*y+20*m.cos(3*t)

def Dx(x,y,t):
    return y


# Alkuarvot ajanhetkellä t=0
y_0 = -8
x_0 = 0.4


# Näytetäänkö tämä kuvaaja? True=Kyllä, False=Ei, huomaa iso alkukirjain!
faasikaavio           = False
funktion_x_kuvaaja    = True
funktion_y_kuvaaja    = False

# Askelkoko, mitä pienempi sitä tarkempi kuvaaja
dt          = 0.0001

# Askelien määrä, mitä suurempi, sitä pidempi väli kuvataan
askelia     = 200000

Numeerinen menetelmä ominaisarvotehtäville

Tässä kappaleessa esitetään eräs menetelmä tavallisen DY:n ominaisarvotehtävän ratkaisemiseksi matriisin ominaisarvojen ja funktion differentiaaliapproksimaation avulla.

Menetelmän idea

Kahdesti derivoituvan funktion ensimmäinen ja toinen derivaatta voidaan laskea myös symmetrisinä raja-arvoina $\begin{align*} y'(x)&=\lim_{h\to0}\frac {y(x+h)-y(x-h)}{2h}\\ y''(x)&=\lim_{h\to0}\frac {y(x+h)-2y(x)+y(x-h)}{h^2}. \end{align*}$ Funktiota y(x) voidaan arvioida kahdella vektorilla ja Vektoriin tulee haluttu väli jaettuna tasaisesti tarkkuuden määräämään palaseen. Vektoriin tulee funktion arvot vastaavissa $X\text:$ n pisteissä. Esimerkiksi funktio x^2 arvioituna välillä [-1,2] tarkkuudella olisi $\begin{align*} X&=(-1,-\tfrac12,0,\tfrac12,1,\tfrac32,2)=(x_1,x_2,...,x_N)\\ Y&=(~~1,~~\tfrac14\,,0,\tfrac14,1,\tfrac94,4)=(y(x_1),y(x_2),...,y(x_N)). \end{align*}$ Voimme laskea vektorin "derivaatan" käyttämällä derivaattojen raja-arvoista saatavia approksimaatioita $\begin{align*} y'(x_n)&\approx\frac1{2h}y(x_{n+1})-\frac1{2h}y(x_{n-1})\\ y''(x_n)&\approx \frac{1}{h^2}(y(x_{n+1})-2y(x_n)+y(x_{n-1})). \end{align*}$ Ensimmäinen derivaatta vastaa kertomista matriisilla $D_1*Y= \frac1{2h} \begin{vmatrix} 0 & 1 & 0 & 0 &...\\ -1 & 0 & 1 & 0 &...\\ 0 &-1 &0 &1 &...\\ 0 &0 &-1& 0 &...\\ \vdots & \vdots &\vdots &\vdots & \end{vmatrix} *Y=Y'.$ Matriisitulosta saatu vektori vastaa funktion derivaattaa kaikkialla paitsi määrittelyvälin päissä.

Edelleen voimme arvioida differentiaalioperaattoria O(y)=ay''+by'+cy matriisitulolla O(Y)=(aD_2+bD_1+cI)*Y=M*Y, missä on identtinen matriisi. Jos haluamme korvata vakiot a,b ja funktioilla $f(x),\,g(x)$ ja h(x), pitää matriisin diagonaaleilla olevat luvut korvata funktioiden arvoilla niitä vastaavissa $X\text:$ n pisteissä.

Ominaisarvotehtävää $O(y)=\lambda y$ vastaa matriisin ominaisarvotehtävä $M*Y=\lambda Y.$ Matriisin ominaisarvoilla voidaan arvioida differentiaalioperaattorin ominaisarvoja ja ominaisvektoreilla niitä vastaavia ominaisfunktioita. Reunaehtona ratkaisufunktioiden on oltava päissä nollaa, mikä seuraa derivaattojen approksimaation "katkeamisesta" määrittelyvälin päissä?

Sovellus: Toteutus Octavella

Seuraavalla ohjelmanpätkällä voi arvioida numeerisesti alkuarvotehtävän $\begin{cases} f(x)\red y''(\green x)+g(x)\red y(\green x)+U(x)\red y(\green x)=\lambda \red y(\green x)\\ \red y(a)=0\\ \red y(b)=0 \end{cases}$ ominaisarvoja $\lambda$ sekä näitä vastaavia ominaisfunktioita.

# dyominaisarvonumer

Värähtelevä kieli numeerisesti

Ratkaistaan ominaisarvotehtävä $\begin{cases} y''(x)=\lambda y\\ y(0)=0\\ y(1)=0 \end{cases}$ numeerisesti. Syöttämällä tiedot ( $a=0,\;b=1,\;f(x)=1$ ) ylläolevaan koodiin saadaan ensimmäiseksi kuudeksi ominaisarvoksi $\begin{align*} -9.8311 && -39.3244 && -88.4795\\ -157.2961 && -245.7734 && -353.9106. \end{align*}$ Alempana tehtävä on laskettu analyyttisesti, jolloin ominaisarvoiksi saadaan $\lambda=-n^2\pi^2.$ Pyöristämällä nämä tarkat arvot neljän desimaalin tarkkuuteen voimme arvioida numeerisen menetelmän antaman tuloksen virhettä. Tarkat arvot ovat pyöristettynä ovat $\begin{align*} -9.8696 && -39.4784 && -88.8264\\ -157.9137 && -246.7401 && -355.3056. \end{align*}$

Kvanttimekaaninen värähtelijä numeerisesti

Schrödingerin yhtälö kvanttimekaaniselle harmoniselle värähtelijälle on $-\frac{\hbar^2}{2m}\red \psi''(\green x)+\frac{1}{2}m\omega^2\green x^2\red \psi(\green x)=\blue E\red \psi(\green x).$ Selvitetään värähtelyn nollapiste-energia, eli pienin ominaisarvo $\blue E.$

Reunaehtoina ratkaisufunktion on mentävä nollaan kun pienenee tai kasvaa rajatta. Numeerista arviointia varten reunaehdoiksi riittää asettaa tarpeeksi suuri luku kuten $\red \psi(-5)=\red \psi(5)=0.$

$\hbar=6$ on redusoitu planckin vakio
on redusoitu massa
on värähtelyn taajuus

Sijoittamalla nämä saadaan $-9\psi''+9x^2\psi=E\psi.$ Numeerinen algoritmi antaa pienimmäksi ominaisarvoksi 8.999. Analyyttinen ratkaisu antaa nollapiste-energialle kaavan $E=\frac12\hbar w=\frac126\cdot3=9.$

tarkkuus = 2^10;
a = -5;
b = 5;

%Jaetaan annettu väli tasaisesti paloihin tarkkuuden mukaan
h=(b-a)/(tarkkuus-1);
x=(a:h:b);

%Määritetään DY:n kertoimet:
%Potentiaalifunktio, eli y:n kerroin
U_x = 9*x.^2;
%Ensimmäisen derivaatan y' kerroin
g_x = 0;
%Toisen derivaatan y'' kerroin
f_x = -9;

Morsepotentiaali ja kaksiatomisen molekyylin värähtely

Morsepotentiaali $V(R)=D_e(1-e^{-\beta(R-R_e)})^2$ Schrödingerin yhtälö $-\frac{\hbar^2}{2m}\psi''+V(R)\psi=E\psi$

De=0.190;
beta=0.9272;
re=2.4086;

tarkkuus = 2^9;
a = re-log(3)/beta;
b = re-log(0.01)/beta;

%Jaetaan annettu väli tasbetasesti paloihin tarkkuuden mukaan
h=(b-a)/(tarkkuus-1);
x=(a:h:b);




%Määritetään DY:n kertoimet:
%Potentiaalifunktio, eli y:n kerroin
U_x = U_x = De*(1-exp(-beta*(x-re))).^2;
%Ensimmäisen derivaatan y' kerroin
g_x = 0;
%Toisen derivaatan y'' kerroin
f_x = -2.7985e-04;

# reuna-arvotehtavat

Reuna-arvotehtävistä

Reuna-arvotehtävässä on DY:n lisäksi annettu funktion $\red y (\green x)$ tai sen jonkin derivaatan arvo useammassa pisteessä.

VLHDY:n ratkaisumetodit on esitelty osassa 1, jonka pohjalta ratkaisut annetaan. Tässä tutkitaan vain, miten annetut reunaehdot vaikuttavat ratkaisun löytymiseen. Reuna-arvot voivat muuttaa annettujen ehtojen toteuttavien funktioiden lisäksi ratkaisujen lukumäärää.

Esimerkki: Reuna-arvotehtävä $\begin{cases} \red y''(\green x)+\red y(\green x)=0\\ \red y(\green 0)=\blue 0\\ \red y(\green{\frac{\pi}{2}})= \blue 2 \end{cases}$

Tämä on VLHDY, joka ratkeaa karakteristisen yhtälön avulla. Ratkaisuiksi saadaan $\red y(\green x)=K_1\sin(\green x)+K_2\cos(\green x).$ Ensimmäisestä reuna-arvosta saadaan $\begin{aligned} \red y(\green 0)&=K_1\sin(\green 0)+K_2\cos(\green 0)\\ &=K_2= \blue 0 \end{aligned}$ ja toisesta $\begin{aligned} \red y\Big(\green {\frac{\pi}{2}}\Big)&=K_1\sin\Big(\green {\frac{\pi}{2}}\Big)+K_2\cos\Big(\green {\frac{\pi}{2}}\Big)\\ &=K_1=\blue 2. \end{aligned}$ Siispä reuna-arvotehtävällä on yksi ratkaisu $\red y(\green x)=2\sin(\green x).$ Erilaisilla reuna-arvoilla saadaan erilaisia (ja eri määrä!) ratkaisuja. Esimerkiksi reunaehdoista $\begin{cases} \red y(\green 0)=\blue 0\\ \red y(\green \pi)=\blue 2 \end{cases}$ seuraa, että K_2=0 ja K_2=-2 joten näillä reuna-arvoilla tehtävällä ei ole ratkaisua! Toisaalta ehdoilla $\begin{cases} \red y(\green 0)=\blue 2\\ \red y(\green \pi)=\blue -\blue 2 \end{cases}$ saamme ratkaisuiksi kaikki funktiot muotoa $\red y(\green x)=K_1\sin(\green x)+2\cos(\green x).$

# ominaisarvotehtavat

Ominaisarvotehtävistä

Reuna-arvotehtävää kutsutaan ominaisarvotehtäväksi silloin, kun

DY on lineaarinen ja homogeeninen $a(\green x)\red y''(\green x)+b(\green x)\red y'(\green x)+c(\green x)\red y(\green x)=0$
ja termillä $\red y(\green x)$ on lisäksi tuntematon vakiokerroin $a(\green x)\red y''(\green x)+b(\green x)\red y'(\green x)+c(\green x)\red y(\green x)=\lambda \red y(\green x)$ $\lambda$ .

Tällöin jokainen ratkaisu $\red y(\green x)$ on DY:n ominaisfunktio, jota vastaa tietty kerroin $\lambda,$ eli ominaisarvo.

Vertaus matriisien ominaisarvoihin

$\lambda$ on matriisin A ominaisarvo ja $\overline v$ sitä vastaava ominaisvektori, jos $A\ve{v}=\lambda \ve v.$ Matriisilla kertominen (tai operointi) ei siis muuta ominaisvektorin $\ve v$ suuntaa vaan ainoastaan pituutta (negatiiviset ominaisarvot kääntävät vektorin).

Matriisin voi ajatella operaattorina, joka muuntelee vektoreita. Määritellään differentiaalioperaattori , joka muuntelee funktioita $\red y(\green x)$ . $O(\red y(\green x))=a(x)\frac{d^2}{dx^2}\red y(\green x)+b(x)\frac{d}{dx}\red y(\green x)+c(x)\red y(\green x).$ Operaattorille vastaava ominaisarvotehtävä on $O(\red y(\green x))=\lambda \red y(\green x).$ Tämä voidaan kirjoittaa differentiaaliyhtälöksi $a(x)\red y''(\green x)+b(x)\red y'(\green x)+c(x)\red y(\green x)=\lambda \red y(\green x).$ Tämän yhtälön ratkaisut ovat operaattorin ominaisfunktiota ja lambdan arvot ominaisarvoja.

Fysiikassa tällaisia ongelmia esiintyy muun muassa värähtelyn ja lämmön johtumisen yhteydessä. Myös kemiassa esimerkiksi Schrödingerin yhtälöstä muotoutuu ominaisarvo-ongelmia DY:iden muodossa.

Esimerkki : Ominaisarvotehtävä $\begin{cases} \red y''(\green x)=\blue \lambda \red y(\green x)\\ \red y(\green 0)=\blue 0\\ \red y(\green l)=\blue 0 \end{cases}$

Lähdetään tarkastelemaan minkälaisia ratkaisuja VLHDY $\red y''(\green x)-\blue \lambda \red y(\green x)=0$ saa eri lambdan arvoilla.

Jos $\blue \lambda>0$ saamme rakaisuiksi $\red y(\green x)=K_1e^{\sqrt{\blue \lambda}\green x}+K_2e^{-\sqrt{\blue \lambda}\green x}.$ Koska $e^{\green x} \neq 0$ kaikilla $\green x$ , on ainoa reunaehdot toteuttava ratkaisu siis , eli $\red y(\green x)=0$ kaikilla $\green x.$
Jos $\blue \lambda=0$ saamme DY:n ratkaisuiksi $\red y(\green x)=K_1\green x+K_2.$ Myös tässä tapauksessa reunaehdot pakottavat funktion nollafunktioksi.
Jos $\blue \lambda<0$ saamme ratkaisuiksi (huomaa: $-\blue \lambda >0$ ) $\red y(\green x)=K_1\cos(\green x\sqrt{-\blue \lambda})+K_2\sin(\green x\sqrt{-\blue \lambda}).$ Reunaehdoista seuraa, että $\begin{cases} \red y(\green 0)=K_1\underbrace{\cos(0)}_{=1}+K_2\underbrace{\sin(0)}_{=0}=0\\ \red y(\green l)=K_1\cos(l\sqrt{-\blue \lambda})+K_2\sin(l\sqrt{-\blue \lambda} )=0, \end{cases}$ eli $\begin{cases} K_1=0\\ K_2\sin(l\sqrt{-\blue \lambda})=0. \end{cases}$ Jos saisimme jälleen nollafunktion, mutta nyt voimme löytää myös muita ratkaisuja sinin nollakohdista. $\sin(l\sqrt{-\blue \lambda})=0\iff l\sqrt{-\blue \lambda}=\blue n\pi$

Ominaisarvotehtävän ratkaisu näillä reunaehdoilla on $\red y_{\blue n}(\green x) = K_2\sin\left(\frac{\blue n\pi}{l}\green x\right)$ ja ominaisarvot ovat $\blue \lambda=-\dfrac{\blue n^2\pi^2}{l^2},$ kun $\blue n=0,1,2,...$ Tämä tehtävä tulee esiintymään monissa eri yhteyksissä fysiikassa, kuten lämpöyhtälössä, hiukkasen aaltoyhtälössä ja värähtelyissä.

# sarjaratkaisut

Sarjaratkaisuista

Lineaariselle differentiaaliyhtälölle $a(\green x)\red y''(\green x)+b(\green x)\red y'(\green x)+c(\green x)\red{y}(\green x)=0$ voidaan etsiä ratkaisua potenssisarjojen avulla. Tämä tarkoittaa sitä, että sijoitetaan yhtälöön ratkaisu muodossa $\red y(\green x)=\sum_{n=0}^\infty a_n\green x^n$ ja katsotaan mitä ehtoja DY asettaa kertoimille.

Jos DY:ssä $a(\green 0)=0,$ tilanne mutkistuu hieman ja ratkaisua kannattaa etsiä muodossa $\red y(\green x)=\green x^s\sum_{n=0}^\infty a_n\green x^n.$

Sarjaratkaisuihin johtavia erikoistapauksia

Nimettyjä differentiaaliyhtälötyyppejä, joille löytyy sarjaratkaisuja, ovat esimerkiksi

Hermiten DY:t $\red y''-2\green x\red y'+2\lambda \red y=0,$ jonka ratkaisut ovat Hermiten polynomeja.
Besselin DY:t $x^2\red y''+\green x\red y'+(\green x^2-a^2)\red y=0,$ jonka ratkaisut ovat Besselin funktioita.
Legendren DY:t $(1-\green x^2)\red y''-2\green x \red y'+n(n+1)\red y=0,$ jonka ratkaisut ovat Legendren polynomeja ja liittopolynomeja.

Esimerkki 1: Sarjaratkaisu DY:lle $y''(x)-\blue{2x}y'(x)+ \red 2 y(x)=0$

Etsitään ratkaisua potenssisarjalla. Yritefunktio ja sen derivaatat ovat

$\begin{align*} y(x)&=\sum_{n=0}^{\infty}\red {a_n}x^n,\\ y'(x)&=\sum_{n=1}^{\infty}\blue{na_n}x^{n-1}=\sum_{n=0}^{\infty}\blue{na_n}x^{n-1},\\ {y''}(x)&=\sum_{n=2}^{\infty}\green{n(n-1)a_n}x^{n-2}\\ &=\sum_{n=0}^{\infty}\green{(n+1)(n+2)a_{n+2}}x^{n}. \end{align*}$

Ensimmäiseen derivaattaan otetaan termi n=0 mukaan helpottamaan sievennystä. Toisessa derivaatassa vaihdetaan summausindeksi alkamaan nollasta. Sijoittamalla nämä DY:hyn saadaan $\begin{align*} 0=\sum_{n=0}^{\infty}&\green{(n+1)(n+2)a_{n+2}}x^{n}\\ & \blue{-2x}\sum_{n=0}^{\infty}\blue{na_n}x^{n-1}+\red 2\sum_{n=0}^{\infty}\red {a_n}x^n\\ =\sum_{n=0}^{\infty}&(\green{(n+1)(n+2)a_{n+2}} \blue{-2na_n}+\red{2a_n})x^n. \end{align*}$ Tämä on nolla kaikilla vain jos jokaisen termin kerroin on nolla, eli $\green{(n+1)(n+2)a_{n+2}} \blue{-2na_n}+\red{2a_n}=0.$ Tästä saamme kertoimille rekursioyhtälön $\green {a_{n+2}}=\frac{\violet 2(\blue n\red{- 1})}{\green{(n+1)(n+2)}}\violet{a_n}.$ Luku a_0 määrittää kaikki parilliset kertoimet $a_{2n}.$

Tällä saadaan ratkaisu $\begin{aligned} y_1(x)&=a_0\sum_{n=0}^{\infty}(0-1)(2-1)...(2n-3)\frac{2^nx^{2n}}{(2n)!}\\ &=a_0\sum_{n=0}^{\infty}\frac{2^nx^{2n}}{(2n)!}\prod_{i=0}^{n-1}(i-1). \end{aligned}$

Toinen lineaarisesti riippumaton ratkaisu saadaan määrittämällä kerroin a_1 . Huomataan, että rekursioyhtälö pakottaa kaikki muut parittomat kertoimet nollaksijos tämä tuntuu epäselvältä yritä laskea a_3 ja a_5 rekursioyhtälöstä. DY:n ratkaisukannan toinen funktio on siis y_2(x)=x. Kaikki ratkaisut saadaan summamalla $y(x)=K_1y_1(x)+K_2y_2(x) =K_2x+K_1\sum_{n=0}^{\infty}\frac{2^nx^{2n}}{(2n)!}\prod_{i=0}^{n-1}(i-1).$

Esimerkki 2: Edellinen ominaisarvo-ongelmana ja Hermiten polynomit $y''(x)\blue{-2x}y'(x)=\red \lambda y(x)$

Millä lambdan arvoilla DY:llä on polynomiratkaisuja?

Valitaan yritteeksi potenssisarja $y(x)=\displaystyle\sum_{n=0}^{\infty} a_nx^{n},$ tämän derivaatat saamme esimerkistä 1.

Sijoittamalla saadaan $\begin{aligned} 0=&\sum_{n=0}^{\infty}\green{(n+1)(n+2)a_{n+2}}x^n\\ &\blue {-2x}\sum_{n=0}^{\infty}\blue{ na_n}x^{n-1} \red{-\lambda}\sum_{n=0}^{\infty}\red{ a_n}x^n.\\ =&\sum_{n=0}^{\infty}(\green{(n+1)(n+2)a_{n+2}}-\blue{2na_n}-\red{\lambda a_n})x^n \end{aligned}$ Tästä saadaan kerrointen rekursioyhtälö $a_{n+2}=\frac{\blue {2n}+\red \lambda}{\green{(n+1)(n+2)}}a_n.$ Jotta saisimme ratkaisuksi polynomin, on oltava $\blue{ 2n}+\red \lambda=0$ jollakin n=0,1,2,..., jolloin summa katkeaa. Siispä DY:llä on polynomiratkaisuja vain kun $\red\lambda=-\blue{2n}.$ Nämä ratkaisut ovat sopivaa kerrointa vaille Hermiten polynomeja.

# kvanttimekaaninenvarahtelija

Sovellus: Kvanttimekaaninen värähtelijä

Schrödingerin yhtälö kvanttimekaaniselle harmoniselle värähtelijälle on $-\frac{\hbar^2}{2m}\red \psi''(\green x)+\frac{1}{2}m\omega^2\green x^2\red \psi(\green x)=\blue E\red \psi(\green x),$ eli $\red \psi''(\green x)+\left(\frac{2m}{\hbar^2}\blue E-\frac{m^2\omega^2}{\hbar^2}\green x^2\right)\red \psi(\green x)=0,$ missä

$\hbar,$ ja $\omega$ ovat vakioita,
$\red\psi(\green x)$ on tuntematon aaltofunktio,
$U(\green x)=\dfrac{1}{2}m\omega^2\green x^2$ on harmoninen potentiaali ja
$\blue E$ on kokonaisenergia.

Haluamme löytää yhtälölle ratkaisuja, jotka voidaan normalisoida, eli vähintään $\red \psi(\green x)\to 0,~~ \text{ kun }~~ \green x \to \pm\infty.$ Tehdään muuttujanvaihto $\green y=\sqrt{\frac{m\omega}{\hbar}}\green x,$ jolloin $\red \psi''(\green y)+(\frac{2\blue E}{\hbar \omega}-\green y^2)\red \psi(\green y)=0.$ Kun $\green y$ on suuri, aaltoyhtälön ratkaisut käyttäytyvät melkein kuin $\red \psi(\green y)\approx K_1\exp\bigg(\frac{-\green y^2}{2}\bigg)+K_2\exp\bigg(\frac{\green y^2}{2}\bigg).$ Funktio on normitettavissa vain jos K_2=0 . Etsitään yhtälön ratkaisuja yritteellä $\red \psi(\green y)=\violet u(\green y)\exp\left(\frac{-\green y^2}{2}\right).$ Derivoidaan tätä kahdesti, sijoitetaan ja sievennetään, jolloin saadaan DY $\violet{u''}(\green y)-2\green y\violet {u'}(\green y)+\left(\frac{2\blue E}{\hbar \omega}-1\right)\violet u(\green y)=0.$ Kun $\dfrac{2\blue E}{\hbar \omega}-1=2\lambda,$ missä $\lambda=0,1,2...,$ saamme ratkaisuksi Hermiten polynomeja $H_\lambda(\green y)$ ; katso Esimerkki 2 yllä. Muilla $\lambda$ :n arvoilla saisimme jotain, mikä ei ole normitettavissa. Värähtelijän normittamattomat aaltofunktiot ovat siis muotoa $\red \psi_\lambda(\green y)=\violet u_\lambda(\green y) \exp\bigg(\dfrac{-\green y^2}{2}\bigg)= H_\lambda\bigg(\sqrt{\frac{m\omega}{\hbar}}\green x\bigg) \exp\bigg(-\dfrac{m\omega}{2\hbar }\green x^2\bigg)$ ja sen mahdolliset energiatilat ovat $\frac{2\blue E}{\hbar \omega}-1=2\lambda \ \ \ \ \ \Leftrightarrow \ \ \ \ \ \blue E=\hbar \omega\left(\lambda+\frac{1}{2}\right).$