AEP 4830 HW4 Ordinary Differential Equations and Chaotic Systems

4th Order Runge Kutta ODE Solver

Suppose we have the a set of \(N\) ODEs with the following form.

\begin{equation} \frac{d\vec{y}}{dt}=\vec{f}(t,\vec{y})\\ \end{equation}

where \(\vec{y}\) and \(\vec{f}\) are N-dimensional vectors of functions. The 4th order Runge Kutta takes the following algorithm each step \(h\) forward and by calculating a large number of steps, we can have the numerical solution to the set of ODEs.(W. H. Press 2007)

\begin{equation} \vec{y_{n}}=\vec{y}(t_{n})\\ \end{equation} \begin{aligned} \vec{k_{1}}=h\vec{f}(t_{n},\vec{y_{n}}) \\ \vec{k_{2}}=h\vec{f}(t_{n}+0.5h,\vec{y_{n}}+0.5\vec{k_{1}}) \\ \vec{k_{3}}=h\vec{f}(t_{n}+0.5h,\vec{y_{n}}+0.5\vec{k_{2}}) \\ \vec{k_{4}}=h\vec{f}(t_{n}+h,\vec{y_{n}}+\vec{k_{3}})\\ \end{aligned}

The value for \(\vec{y_{n+1}}\) one step away can be calculated.

\begin{equation} \vec{y_{n+1}}=\vec{y_{n}}+\frac{1}{6}\left(\vec{k_{1}}+2\vec{k_{2}}+2\vec{k_{3}}+\vec{k_{4}}\right)\\ \end{equation}

Finally, by starting from \(\vec{y_{n+1}}\), the solution can be obtained.

Testing

After constructing the 4th order of Runge Kutta method to solve any general ODEs, it’s safe if we test it on an ODE that has the known solution. The testing function is a simple harmonic oscillator without damping.

\begin{equation} y^{\prime\prime}(t)+\omega^{2}y(t)=0\\ \end{equation}

with initial conditions \(y(0)=1\) and \(y^{\prime}(0)=0\). The solution is obvious \(y(t)=\cos(t)\) In order to implement 4th order Runge Kutta, we can rearrange the ODE into a set of 2 first order ODEs.

\begin{equation} \frac{dy_{0}}{dt}=y_{1};\frac{dy_{1}}{dt}=-\omega^{2}y_{0}\\ \end{equation}

The data were generated into a .dat file and were loaded and plotted using MATLAB.
The testing solution is plotted in Fig. (1) with amplitude and period quantitatively printed in the .exe window, as shown in Fig. (2). In addition, the phase space plot, \(y^{\prime}(t)\) vs \(y(t)\) is shown in Fig. (3). The period is found to be \(T=3.13\) which is supposed to be \(\pi\) and The amplitude is then \(A=0.999965\). The errors here are because of the finite step size \(h=0.002\). The phase space plot is a circle corresponding to the following parametric functions.

\begin{equation} y=\sin(t);x=\cos(t)\\ \end{equation}

The solution is stable and periodic, which verifies the feasibility of the 4th order Runge Kutta code.

The numerical solution of simple harmonic oscillation using 4th order Runge Kutta code. The step size is \(h=0.002\) and the number of steps is \(nstep=10000\).