Differentialekvationer kan lösas numeriskt i Microsoft Excel med hjälp av exempelvis Eulers stegmetod:

y_{n+1}=y_n+y_n^\prime\ast\Delta t

Första ordningens ordinära differentialekvationer är på formen

y^\prime+ay=0,

medan andra ordningens ordinära differentialekvationer är på formen

y^{\prime\prime}+ay^\prime+by=0

vilket innebär att vid beräkning av andra ordningens ordinära differentialekvationer med Eulers stegmetod behövs även beräkningssteget

y_{n+1}^\prime=y_n^\prime+y_n^{\prime\prime}\ast\Delta t

Här kommer att lösas de tre rörelseekvationerna som beskriver svängningar för en enkel mekanisk modell av koldioxidmolekylen, där Langrangian och Euler-Lagranges ekvationer ger följande rörelseekvationer för de tre i koldioxiden ingående atomerna:

x_1:\ \ \ \ m{\ddot{x}}_1-k\left(x_2-x_1\right)=0 ,

x_2:\ \ \ \ M{\ddot{x}}_2-k\left(x_1-2x_2+x_3\right)=0 ,

x_3:\ \ \ \ m{\ddot{x}}_3-k\left(x_2-x_3\right)=0\

För att kunna använda Eulers stegmetod skrivs dessa om enligt:

{\ddot{x}}_1=\frac{k}{m}\left(x_2-x_1\right) ,

{\ddot{x}}_2=\frac{k}{M}\left(x_1{-2x}_2+x_3\right) ,

{\ddot{x}}_3=\frac{k}{m}\left(x_2-x_3\right)

Om dessa jämförs med en generell andra ordningens differentialekvation x^{\prime\prime}+ax^\prime+bx=0 framgår att termen x^\prime=\dot{x} saknas i de tre rörelseekvationerna, men termen behövs för att kunna använda Eulers stegmetod för andra ordningens differentialekvationer.

Vid konstant acceleration gäller dock att \dot{x}=v=at där a=\ddot{x}\ \ \rightarrow\ \ \ \dot{x}=\ddot{x}\ast\Delta t ,
och på det korta intervallet \Delta t som kommer att användas i Eulers stegmetod kan accelerationen anses vara konstant.

För differentialekvationer behövs även begynnelsevillkor och för den mekaniska modellen utgörs dessa av:

x_1,\ x_2 och x_3.

Begynnelsevillkoren sätter den mekaniska modellen i svängning och lämpliga begynnelsevärden kan inte sättas hur som helst. Den verkliga koldioxidmolekylen svänger harmoniskt enligt de egenfrekvenser som hör ihop med det vågtal som beskriver en viss svängsningmode, så därav bör även den mekaniska modellen sättas i harmonisk svängning. Genom att utgå ifrån storleken på koldioxidmolekylen kan lämpliga begynnelsevillkor för harmonisk svängning hittas.

Om endast kolatomen flyttas något åt ettdera hållet bör harmonisk svängning uppstå och en lämplig distans kan vara 50 pm.

När kolatomen flyttas åt höger kommer fjädern till vänster om kolatomen att sträckas ut och fjädern till höger om kolatomen att tryckas ihop. När svängningen kommer igång från begynnelsevärdena,

x_1=0,\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ x_2=5\ast{10}^{-11}\ m,\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ x_3=0,

kommer fjädrarna att flytta kolatomen till vänster, vänstra syreatomen till höger (vänstra fjädern drar ihop sig) och högra syreatomen till höger (högra fjädern trycker ut sig). Genom att beräkna de tre ordinära differentialekvationerna fås reda på hur dessa svängningar kommer att gestalta sig. Innan genomgång av Excelberäkningar sker visas resultatet av beräkningarna, där det tydligt framgår att valda begynnelsevillkor innebär att alla tre atomerna i koldioxidmolekylen svänger harmoniskt:

Ur graferna kan utläsas att från begynnelsevillkoren (startvärden vid 1.) kommer kolatomen att svänga tillbaka åt vänster, över viloläget med -22,7 pm, samtidigt som båda syreatomerna svänger +27,3 pm åt höger (2.) Exakta värden är dock avlästa i Exceltabellerna bakom graferna och inte ur graferna direkt. Sedan (3.), svänger atomerna tillbaka till tstartläget varvid den cykliska svängningsrörelsen startar om igen och fortsätter på samma sätt.

Ur grafen går även utläsa att periodtiden är 0,124 ps, vilket innebär att det handlar om snabba svängingar (vilket det alltid handlar om på atomnivå).

Denna periodtiden kan kontrolleras mot den periodtid som fås ur de egenvinkelfrekvenser som svängningarna baseras på enligt:

\omega_2=5,15287\ast{10}^{14}\ rad/s \omega=\frac{2\pi}{T}\ \ \rightarrow\ \ \ T=\frac{2\pi}{\omega_2}=\frac{2\pi}{5,15287\ast{10}^{14}}=1,22\ast{10}^{-14}=0,122\ ps

Skillnaden mellan 0,122 ps och 0,124 ps rör sig om 1,6 procent och kan delvis bero på något avrundade siffror i de numeriska beräkningarna. Men numerisk beräkning av differentialekvationer vid snabba svängningar är kritiska och mer avancerade programvaror än Excel, har lite olika metoder att handskas med detta problem. I Excel är det möjligt att minimera problemet, vilket beskrivs under genomgången av Excelberäkningarna av differentialekvationerna.

Skapa ett Excelark med de kolumner och rader som behövs.

Delta t (dt) kan ses som ”upplösningen” för beräkningarna, det vill säga hur ofta beräkning ska ske. Här väljs dt drygt en faktor tusen mindre än periodtiden vilket innebär att under en period kommer 1220 beräkningar att ske, som innebär lika många punkter på linjen/sinuskurvan under. Om ”upplösningen” är för låg blir beräkningsresultatet felaktigt och beräkningen av differentialekvationerna skenar fort iväg åt något håll.

\frac{T}{dt}=\frac{1,22\ast{10}^{-14}}{1\ast{10}^{-17}}=1220

Första kolumnen för iterationer är inte helt nödvändig, men det är ett bra sätt att se hur många iterationer eller beräkningssteg som görs. Eftersom första raden egentligen inte är en iteration kan väljas att ge den iterationsnumret 0, vilket även innebär att numreringen ”följer” tiden numrärt. För under kolumnen tid adderas från 0, dt (delta t) för varje rad (eller iteration), där det finns ett flertal sätt i Excel att göra detta. Ett sätt är att använda värdet från raden ovan och addera till dt, där cellen M1 innehållande dt måste låsas. Genom att använda sig av detta sätt är det enkelt att ändra och laborera med olika värden på dt.

Ifyllnad av första raden

Till cellerna för x1, x2 och x3 tas värden från cellerna O1, O2 och O3, där dessa celler måste låsas. Sedan beräknas övriga celler på raden enligt nedan, där de värden som tas från celler i kolumn M måste låsas.

För andra raden används Eulers stegmetod
y_{n+1}=y_n+y_n^\prime\ast\Delta t,\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ y_{n+1}^\prime=y_n^\prime+y_n^{\prime\prime}\ast\Delta t

Till skillnad från {\ddot{x}}_1, {\ddot{x}}_2 och {\ddot{x}}_3, som beräknas utifrån värden på raden tillsammans med låsta värden från M-kolumnen, så använder Eulers stegmetod värden från ovanstående rad. Vilket innebär följande beräkningar:

{\ddot{x}}_1, {\ddot{x}}_2 och {\ddot{x}}_3, beräknas på samma sätt som raden ovanför, varvid exempelvis funktionen ”fyll nedåt” kan användas, alternativt kopieras ”formeln” (ej bara värdet) från raden ovan och klistras in.

Tredje raden

Tredje raden och resterande rader kan sedan beräknas likadant som andra raden, genom att markera samtliga kolumner på andra raden (nummer 3 i Excelarket) och fylla nedåt. Fyll nedåt görs genom att fatta tag i kvadraten längs ner till höger vid de markerade rederna och dra den nedåt de önskade antalet rader. Eftersom en period rör sig om 1220 rader kan ett par tusen rader fyllas för att få med flera perioder.

Differentialekvationerna är nu beräknande och plottar på de tre x-kolumnerna som funktion av tiden kan skapas för att se resultatet. Jag väljer dock istället att beräkna resterande rader genom att utöka Eulers enstegsmetod till en tvåstegsmetod, vilket ökar precisionen i beräkningarna. Eulers enstegsmetod använder värden från den senaste beräkningen (raden ovanför/det tidigare steget) för att beräkna aktuellt steg/rad.
Adams-Bashforth är en tvåstegsmetod som använder värden från de två senaste beräkningarna (från två tidigare steg) för att beräkna aktuellt steg/rad enligt:

x_{n+2}=x_{n+1}+{1,5x}_{n+1}^\prime\ast\Delta t-{0,5x}_n^\prime\ast\Delta t\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ x_{n+2}^\prime=x_{n+1}^\prime+1,5x_{n+1}^{\prime\prime}\ast\Delta t-0,5x_n^{\prime\prime}\ast\Delta t
{\ddot{x}}_1, {\ddot{x}}_2 och {\ddot{x}}_3 beräknas på samma sätt som raden ovanför och avslutningsvis markeras samtliga kolumner på tredje raden (nummer 4 i Excelarket) och funktionen fyll nedåt körs i några tusen rader, här är 15 002 rader valda.

Plottar

x_1,\ x_2\ och x_3 kan nu plottas som funktion av tiden genom att markera de kolumner som ska plottas, samt därefter välja infoga ett lämpligt diagram.

Dessa differentialekvationer är beräknade för en koldioxidmolekyl som svänger med vågtalet,

\bar{v}=2734\ cm^{-1}, som hör ihop med egenvinkelfrekvensen \omega=5,15287\ast{10}^{14}\ rad/s .

Läs mer under beräkna egenvinkelfrekvenser för koldioxidmolekylen.