The order of the ODE refers to the order of the derivative on the left-hand side in the equation $$ \begin{equation} \frac{dy}{dt}=f(t,y). \label{_auto1} \end{equation} $$ This equation is of first order and \( f \) is an arbitrary function. A second-order equation goes typically like $$ \begin{equation} \frac{d^2y}{dt^2}=f(t,\frac{dy}{dt},y). \label{_auto2} \end{equation} $$ A well-known second-order equation is Newton's second law $$ \begin{equation} m\frac{d^2x}{dt^2}=-kx, \label{eq:newton} \end{equation} $$ where \( k \) is the force constant. ODE depend only on one variable
partial differential equations like the time-dependent Schr\"odinger equation $$ \begin{equation} i\hbar\frac{\partial \psi({\bf x},t)}{\partial t}= -\frac{\hbar^2}{2m}\left( \frac{\partial^2 \psi({\bf r},t)}{\partial x^2} + \frac{\partial^2 \psi({\bf r},t)}{\partial y^2}+ \frac{\partial^2 \psi({\bf r},t)}{\partial z^2}\right) + V({\bf x})\psi({\bf x},t), \label{_auto3} \end{equation} $$ may depend on several variables. In certain cases, like the above equation, the wave function can be factorized in functions of the separate variables, so that the Schroedinger equation can be rewritten in terms of sets of ordinary differential equations. These equations are discussed in chapter 10. Involve boundary conditions in addition to initial conditions.
We distinguish also between linear and non-linear differential equation where for example $$ \begin{equation} \frac{dy}{dt}=g^3(t)y(t), \label{_auto4} \end{equation} $$ is an example of a linear equation, while $$ \begin{equation} \frac{dy}{dt}=g^3(t)y(t)-g(t)y^2(t), \label{_auto5} \end{equation} $$ is a non-linear ODE.
Another concept which dictates the numerical method chosen for solving an ODE, is that of initial and boundary conditions. To give an example, if we study white dwarf stars or neutron stars we will need to solve two coupled first-order differential equations, one for the total mass \( m \) and one for the pressure \( P \) as functions of \( \rho \) $$ \frac{dm}{dr}=4\pi r^{2}\rho (r)/c^2, $$ and $$ \frac{dP}{dr}=-\frac{Gm(r)}{r^{2}}\rho (r)/c^2. $$ where \( \rho \) is the mass-energy density. The initial conditions are dictated by the mass being zero at the center of the star, i.e., when \( r=0 \), yielding \( m(r=0)=0 \). The other condition is that the pressure vanishes at the surface of the star.
In the solution of the Schroedinger equation for a particle in a potential, we may need to apply boundary conditions as well, such as demanding continuity of the wave function and its derivative.
In many cases it is possible to rewrite a second-order differential equation in terms of two first-order differential equations. Consider again the case of Newton's second law in Eq. \eqref{eq:newton}. If we define the position \( x(t)=y^{(1)}(t) \) and the velocity \( v(t)=y^{(2)}(t) \) as its derivative $$ \begin{equation} \frac{dy^{(1)}(t)}{dt}=\frac{dx(t)}{dt}=y^{(2)}(t), \label{_auto6} \end{equation} $$ we can rewrite Newton's second law as two coupled first-order differential equations $$ \begin{equation} m\frac{dy^{(2)}(t)}{dt}=-kx(t)=-ky^{(1)}(t), \label{eq:n1} \end{equation} $$ and $$ \begin{equation} \frac{dy^{(1)}(t)}{dt}=y^{(2)}(t). \label{eq:n2} \end{equation} $$
These methods fall under the general class of one-step methods. The algoritm is rather simple. Suppose we have an initial value for the function \( y(t) \) given by $$ \begin{equation} y_0=y(t=t_0). \label{_auto7} \end{equation} $$ We are interested in solving a differential equation in a region in space \( [a,b] \). We define a step \( h \) by splitting the interval in \( N \) sub intervals, so that we have $$ \begin{equation} h=\frac{b-a}{N}. \label{_auto8} \end{equation} $$ With this step and the derivative of \( y \) we can construct the next value of the function \( y \) at $$ \begin{equation} y_1=y(t_1=t_0+h), \label{_auto9} \end{equation} $$ and so forth.
If the function is rather well-behaved in the domain \( [a,b] \), we can use a fixed step size. If not, adaptive steps may be needed. Here we concentrate on fixed-step methods only. Let us try to generalize the above procedure by writing the step \( y_{i+1} \) in terms of the previous step \( y_i \) $$ \begin{equation} y_{i+1}=y(t=t_i+h)=y(t_i) + h\Delta(t_i,y_i(t_i)) + O(h^{p+1}), \label{_auto10} \end{equation} $$ where \( O(h^{p+1}) \) represents the truncation error. To determine \( \Delta \), we Taylor expand our function \( y \) $$ \begin{equation} y_{i+1}=y(t=t_i+h)=y(t_i) + h(y'(t_i)+\dots +y^{(p)}(t_i)\frac{h^{p-1}}{p!}) + O(h^{p+1}), \label{eq:taylor} \end{equation} $$ where we will associate the derivatives in the parenthesis with $$ \begin{equation} \Delta(t_i,y_i(t_i))=(y'(t_i)+\dots +y^{(p)}(t_i)\frac{h^{p-1}}{p!}). \label{eq:delta} \end{equation} $$
We define $$ \begin{equation} y'(t_i)=f(t_i,y_i) \label{_auto11} \end{equation} $$ and if we truncate \( \Delta \) at the first derivative, we have $$ \begin{equation} y_{i+1}=y(t_i) + hf(t_i,y_i) + O(h^2), \label{eq:euler} \end{equation} $$ which when complemented with \( t_{i+1}=t_i+h \) forms the algorithm for the well-known Euler method. Note that at every step we make an approximation error of the order of \( O(h^2) \), however the total error is the sum over all steps \( N=(b-a)/h \), yielding thus a global error which goes like \( NO(h^2)\approx O(h) \).
To make Euler's method more precise we can obviously decrease \( h \) (increase \( N \)). However, if we are computing the derivative \( f \) numerically by for example the two-steps formula $$ f'_{2c}(x)= \frac{f(x+h)-f(x)}{h}+O(h), $$ we can enter into roundoff error problems when we subtract two almost equal numbers \( f(x+h)-f(x)\approx 0 \). Euler's method is not recommended for precision calculation, although it is handy to use in order to get a first view on how a solution may look like. As an example, consider Newton's equation rewritten in Eqs. \eqref{eq:n1} and \eqref{eq:n2}. We define \( y_0=y^{(1)}(t=0) \) an \( v_0=y^{(2)}(t=0) \). The first steps in Newton's equations are then $$ \begin{equation} y^{(1)}_1=y_0+hv_0+O(h^2) \label{_auto12} \end{equation} $$ and $$ \begin{equation} y^{(2)}_1=v_0-hy_0k/m+O(h^2). \label{_auto13} \end{equation} $$
The Euler method is asymmetric in time, since it uses information about the derivative at the beginning of the time interval. This means that we evaluate the position at \( y^{(1)}_1 \) using the velocity at \( y^{(2)}_0=v_0 \). A simple variation is to determine \( y^{(1)}_{n+1} \) using the velocity at \( y^{(2)}_{n+1} \), that is (in a slightly more generalized form) $$ \begin{equation} y^{(1)}_{n+1}=y^{(1)}_{n}+h y^{(2)}_{n+1}+O(h^2) \label{_auto14} \end{equation} $$ and $$ \begin{equation} y^{(2)}_{n+1}=y^{(2)}_{n}+h a_{n}+O(h^2). \label{_auto15} \end{equation} $$ The acceleration \( a_n \) is a function of \( a_n(y^{(1)}_{n}, y^{(2)}_{n},t) \) and needs to be evaluated as well. This is the Euler-Cromer method.
Let us then include the second derivative in our Taylor expansion. We have then $$ \begin{equation} \Delta(t_i,y_i(t_i))=f(t_i)+\frac{h}{2}\frac{df(t_i,y_i)}{dt}+O(h^3). \label{_auto16} \end{equation} $$ The second derivative can be rewritten as $$ \begin{equation} y''=f'=\frac{df}{dt}=\frac{\partial f}{\partial t}+\frac{\partial f}{\partial y}\frac{\partial y}{\partial t}=\frac{\partial f}{\partial t}+\frac{\partial f}{\partial y}f \label{eq:derivatives} \end{equation} $$ and we can rewrite Eq.\ \eqref{eq:taylor} as $$ \begin{equation} y_{i+1}=y(t=t_i+h)=y(t_i) +hf(t_i)+ \frac{h^2}{2}\left(\frac{\partial f}{\partial t}+\frac{\partial f}{\partial y}f\right) + O(h^{3 }), \label{_auto17} \end{equation} $$ which has a local approximation error \( O(h^{3 }) \) and a global error \( O(h^{2}) \).
These approximations can be generalized by using the derivative \( f \) to arbitrary order so that we have $$ \begin{equation} y_{i+1}=y(t=t_i+h)=y(t_i) + h(f(t_i,y_i)+\dots f^{(p-1)}(t_i,y_i) \frac{h^{p-1}}{p!}) + O(h^{p+1}). \label{_auto18} \end{equation} $$ These methods, based on higher-order derivatives, are in general not used in numerical computation, since they rely on evaluating derivatives several times. Unless one has analytical expressions for these, the risk of roundoff errors is large.
The most obvious improvements to Euler's and Euler-Cromer's algorithms, avoiding in addition the need for computing a second derivative, is the so-called midpoint method. We have then $$ \begin{equation} y^{(1)}_{n+1}=y^{(1)}_{n}+\frac{h}{2}\left(y^{(2)}_{n+1}+y^{(2)}_{n}\right)+O(h^2) \label{_auto19} \end{equation} $$ and $$ \begin{equation} y^{(2)}_{n+1}=y^{(2)}_{n}+h a_{n}+O(h^2), \label{_auto20} \end{equation} $$ yielding $$ \begin{equation} y^{(1)}_{n+1}=y^{(1)}_{n}+hy^{(2)}_{n}+\frac{h^2}{2}a_n+O(h^3) \label{_auto21} \end{equation} $$ implying that the local truncation error in the position is now \( O(h^3) \), whereas Euler's or Euler-Cromer's methods have a local error of \( O(h^2) \).
Thus, the midpoint method yields a global error with second-order accuracy for the position and first-order accuracy for the velocity. However, although these methods yield exact results for constant accelerations, the error increases in general with each time step.
One method that avoids this is the so-called half-step method. Here we define $$ \begin{equation} y^{(2)}_{n+1/2}=y^{(2)}_{n-1/2}+h a_{n}+O(h^2), \label{_auto22} \end{equation} $$ and $$ \begin{equation} y^{(1)}_{n+1}=y^{(1)}_{n}+hy^{(2)}_{n+1/2} +O(h^2). \label{_auto23} \end{equation} $$ Note that this method needs the calculation of \( y^{(2)}_{1/2} \). This is done using e.g., Euler's method $$ \begin{equation} y^{(2)}_{1/2}=y^{(2)}_{0}+h a_{0}+O(h^2). \label{_auto24} \end{equation} $$ As this method is numerically stable, it is often used instead of Euler's method.
Another method which one may encounter is the Euler-Richardson method with $$ \begin{equation} y^{(2)}_{n+1}=y^{(2)}_{n}+h a_{n+1/2}+O(h^2), \label{eq:er1} \end{equation} $$ and $$ \begin{equation} \label{eq:er2} y^{(1)}_{n+1}=y^{(1)}_{n}+hy^{(2)}_{n+1/2} +O(h^2). \end{equation} $$
Another set of popular algorithms, which are both numerically stable and easy to implement are the Verlet algorithms, with the velocity Verlet method as widely used in for example Molecular dynamics calculations. Consider again a second-order differential equation like Newton's second law, whose one-dimensional version reads $$ m\frac{d^2 x}{dt^2}= F(x,t), $$ which we rewrite in terms of two coupled differential equations $$ \frac{dx}{dt}=v(x,t) \hspace{1cm}\mathrm{and}\hspace{1cm} \frac{dv}{dt}=F(x,t)/m=a(x,t). $$
If we now perform a Taylor expansion $$ x(t+h) = x(t)+hx^{(1)}(t)+\frac{h^2}{2}x^{(2)}(t)+O(h^3). $$ In our case the second derivative is known via Newton's second law, namely \( x^{(2)}(t)=a(x,t) \). If we add to the above equation the corresponding Taylor expansion for \( x(t-h) \), we obtain, using the discretized expressions $$ x(t_i\pm h) = x_{i\pm 1} \hspace{1cm}\mathrm{and}\hspace{1cm} x_i=x(t_i), $$ we arrive at $$ x_{i+1} = 2x_i - x_{i-1} +h^2x^{(2)}_i+O(h^4). $$
We note that the truncation error goes like \( O(h^4) \) since all the odd terms cancel when we add the two Taylor expansions. We see also that the velocity is not directly included in the equation since the function \( x^{(2)}=a(x,t) \) is supposed to be known. If we need the velocity however, we can compute it using the well-known formula $$ x^{(1)}_i=\frac{x_{i+1}-x_{i-1}}{2h}+O(h^2). $$ We note that the velocity has a truncation error which goes like \( O(h^2) \). In for example so-called Molecular dynamics calculations, since the acceleration is normally known via Newton's second law, there is seldomly a need for computing the velocity.
We note also that the algorithm for the position is not self-starting since, for \( i=0 \) it depends on the value of \( x \) at the fictitious value \( x_{-1} \).
We can amend this by introducing the velocity Verlet method.
We have the Taylor expansion for the position given by $$ x_{i+1} = x_i+hx^{(1)}_i+\frac{h^2}{2}x^{(2)}_i+O(h^3). $$ The corresponding expansion for the velocity is $$ v_{i+1} = v_i+hv^{(1)}_i+\frac{h^2}{2}v^{(2)}_i+O(h^3). $$ Via Newton's second law we have normally an analytical expression for the derivative of the velocity, namely $$ v^{(1)}_i = \frac{d^2 x}{dt^2}\vert_{i}= \frac{F(x_i,t_i)}{m}, $$
If we add to this the corresponding expansion for the derivative of the velocity $$ v^{(1)}_{i+1} = v^{(1)}_i+hv^{(2)}_i+O(h^2), $$ and retain only terms up to the second derivative of the velocity since our error goes as \( O(h^3) \), we have $$ hv^{(2)}_i\approx v^{(1)}_{i+1}-v^{(1)}_i. $$ We can then rewrite the Taylor expansion for the velocity as $$ v_{i+1} = v_i+\frac{h}{2}\left( v^{(1)}_{i+1}+v^{(1)}_{i}\right)+O(h^3). $$
Our final equations for the position and the velocity become then $$ x_{i+1} = x_i+hv_i+\frac{h^2}{2}v^{(1)}_{i}+O(h^3), $$ and $$ v_{i+1} = v_i+\frac{h}{2}\left( v^{(1)}_{i+1}+v^{(1)}_{i}\right)+O(h^3). $$ Note well that the term \( v^{(1)}_{i+1} \) depends on the position at \( x_{i+1} \). This means that you need to calculate the position at the updated time \( t_{i+1} \) before the computing the next velocity. Note also that the derivative of the velocity at the time \( t_i \) used in the updating of the position can be reused in the calculation of the velocity update as well.
Runge-Kutta (RK) methods are based on Taylor expansion formulae, but yield in general better algorithms for solutions of an ODE. The basic philosophy is that it provides an intermediate step in the computation of \( y_{i+1} \).
To see this, consider first the following definitions $$ \begin{equation} \frac{dy}{dt}=f(t,y), \label{_auto25} \end{equation} $$ and $$ \begin{equation} y(t)=\int f(t,y) dt, \label{_auto26} \end{equation} $$ and $$ \begin{equation} y_{i+1}=y_i+ \int_{t_i}^{t_{i+1}} f(t,y) dt. \label{_auto27} \end{equation} $$
To demonstrate the philosophy behind RK methods, let us consider the second-order RK method, RK2. The first approximation consists in Taylor expanding \( f(t,y) \) around the center of the integration interval \( t_i \) to \( t_{i+1} \), that is, at \( t_i+h/2 \), \( h \) being the step. Using the midpoint formula for an integral, defining \( y(t_i+h/2) = y_{i+1/2} \) and \( t_i+h/2 = t_{i+1/2} \), we obtain $$ \begin{equation} \int_{t_i}^{t_{i+1}} f(t,y) dt \approx hf(t_{i+1/2},y_{i+1/2}) +O(h^3). \label{_auto28} \end{equation} $$ This means in turn that we have $$ \begin{equation} y_{i+1}=y_i + hf(t_{i+1/2},y_{i+1/2}) +O(h^3). \label{_auto29} \end{equation} $$
However, we do not know the value of \( y_{i+1/2} \). Here comes thus the next approximation, namely, we use Euler's method to approximate \( y_{i+1/2} \). We have then $$ \begin{equation} y_{(i+1/2)}=y_i + \frac{h}{2}\frac{dy}{dt} = y(t_i) + \frac{h}{2}f(t_i,y_i). \label{_auto30} \end{equation} $$ This means that we can define the following algorithm for the second-order Runge-Kutta method, RK2. $$ \begin{equation} k_1=hf(t_i,y_i), \label{_auto31} \end{equation} $$ $$ \begin{equation} k_2=hf(t_{i+1/2},y_i+k_1/2), \label{_auto32} \end{equation} $$ with the final value $$ \begin{equation} y_{i+i}\approx y_i + k_2 +O(h^3). \label{_auto33} \end{equation} $$
The difference between the previous one-step methods is that we now need an intermediate step in our evaluation, namely \( t_i+h/2 = t_{(i+1/2)} \) where we evaluate the derivative \( f \). This involves more operations, but the gain is a better stability in the solution.
The fourth-order Runge-Kutta, RK4, which we will employ in the solution of various differential equations below, has the following algorithm $$ k_1=hf(t_i,y_i) \hspace{0.5cm} k_2=hf(t_i+h/2,y_i+k_1/2) $$ $$ k_3=hf(t_i+h/2,y_i+k_2/2)\hspace{0.5cm} k_4=hf(t_i+h,y_i+k_3) $$ with the final result $$ y_{i+1}=y_i +\frac{1}{6}\left( k_1 +2k_2+2k_3+k_4\right). $$ Thus, the algorithm consists in first calculating \( k_1 \) with \( t_i \), \( y_1 \) and \( f \) as inputs. Thereafter, we increase the step size by \( h/2 \) and calculate \( k_2 \), then \( k_3 \) and finally \( k_4 \). The global error goes as \( O(h^4) \).
Our first example is the classical case of simple harmonic oscillations, namely a block sliding on a horizontal frictionless surface. The block is tied to a wall with a spring. If the spring is not compressed or stretched too far, the force on the block at a given position \( x \) is $$ F=-kx. $$ The negative sign means that the force acts to restore the object to an equilibrium position. Newton's equation of motion for this idealized system is then $$ m\frac{d^2x}{dt^2}=-kx, $$ or we could rephrase it as $$ \frac{d^2x}{dt^2}=-\frac{k}{m}x=-\omega_0^2x, \label{eq:newton1} $$ with the angular frequency \( \omega_0^2=k/m \).
The above differential equation has the advantage that it can be solved analytically with solutions on the form $$ x(t)=Acos(\omega_0t+\nu), $$ where \( A \) is the amplitude and \( \nu \) the phase constant. This provides in turn an important test for the numerical solution and the development of a program for more complicated cases which cannot be solved analytically.
With the position \( x(t) \) and the velocity \( v(t)=dx/dt \) we can reformulate Newton's equation in the following way $$ \frac{dx(t)}{dt}=v(t), $$ and $$ \frac{dv(t)}{dt}=-\omega_0^2x(t). $$ We are now going to solve these equations using the Runge-Kutta method to fourth order discussed previously.
Before proceeding however, it is important to note that in addition to the exact solution, we have at least two further tests which can be used to check our solution.
Since functions like \( cos \) are periodic with a period \( 2\pi \), then the solution \( x(t) \) has also to be periodic. This means that $$ x(t+T)=x(t), $$ with \( T \) the period defined as $$ T=\frac{2\pi}{\omega_0}=\frac{2\pi}{\sqrt{k/m}}. $$ Observe that \( T \) depends only on \( k/m \) and not on the amplitude of the solution.
In addition to the periodicity test, the total energy has also to be conserved.
Suppose we choose the initial conditions $$ x(t=0)=1\hspace{0.1cm} \mathrm{m}\hspace{1cm} v(t=0)=0\hspace{0.1cm}\mathrm{m/s}, $$ meaning that block is at rest at \( t=0 \) but with a potential energy $$ E_0=\frac{1}{2}kx(t=0)^2=\frac{1}{2}k. $$ The total energy at any time \( t \) has however to be conserved, meaning that our solution has to fulfil the condition $$ E_0=\frac{1}{2}kx(t)^2+\frac{1}{2}mv(t)^2. $$
An algorithm which implements these equations is included below.
The following python program performs essentially the same calculations as the previous c++ code.
The angular equation of motion of the pendulum is given by Newton's equation and with no external force it reads $$ \begin{equation} ml\frac{d^2\theta}{dt^2}+mgsin(\theta)=0, \label{_auto34} \end{equation} $$ with an angular velocity and acceleration given by $$ \begin{equation} v=l\frac{d\theta}{dt}, \label{_auto35} \end{equation} $$ and $$ \begin{equation} a=l\frac{d^2\theta}{dt^2}. \label{_auto36} \end{equation} $$
We do however expect that the motion will gradually come to an end due a viscous drag torque acting on the pendulum. In the presence of the drag, the above equation becomes $$ \begin{equation} ml\frac{d^2\theta}{dt^2}+\nu\frac{d\theta}{dt} +mgsin(\theta)=0, \label{eq:pend1} \end{equation} $$ where \( \nu \) is now a positive constant parameterizing the viscosity of the medium in question. In order to maintain the motion against viscosity, it is necessary to add some external driving force. We choose here a periodic driving force. The last equation becomes then $$ \begin{equation} ml\frac{d^2\theta}{dt^2}+\nu\frac{d\theta}{dt} +mgsin(\theta)=Asin(\omega t), \label{eq:pend2} \end{equation} $$ with \( A \) and \( \omega \) two constants representing the amplitude and the angular frequency respectively. The latter is called the driving frequency.
We define $$ \omega_0=\sqrt{g/l}, $$ the so-called natural frequency and the new dimensionless quantities $$ \hat{t}=\omega_0t, $$ with the dimensionless driving frequency $$ \hat{\omega}=\frac{\omega}{\omega_0}, $$ and introducing the quantity \( Q \), called the quality factor, $$ Q=\frac{mg}{\omega_0\nu}, $$ and the dimensionless amplitude $$ \hat{A}=\frac{A}{mg} $$
We have $$ \frac{d^2\theta}{d\hat{t}^2}+\frac{1}{Q}\frac{d\theta}{d\hat{t}} +sin(\theta)=\hat{A}cos(\hat{\omega}\hat{t}). $$ This equation can in turn be recast in terms of two coupled first-order differential equations as follows $$ \frac{d\theta}{d\hat{t}}=\hat{v}, $$ and $$ \frac{d\hat{v}}{d\hat{t}}=-\frac{\hat{v}}{Q}-sin(\theta)+\hat{A}cos(\hat{\omega}\hat{t}). $$ These are the equations to be solved. The factor \( Q \) represents the number of oscillations of the undriven system that must occur before its energy is significantly reduced due to the viscous drag. The amplitude \( \hat{A} \) is measured in units of the maximum possible gravitational torque while \( \hat{\omega} \) is the angular frequency of the external torque measured in units of the pendulum's natural frequency.
We start with a simpler case first, the Earth-Sun system in two dimensions only. The gravitational force \( F_G \) is $$ F=\frac{GM_{\odot}M_E}{r^2}, $$ where \( G \) is the gravitational constant, $$ M_E=6\times 10^{24}\mathrm{Kg}, $$ the mass of Earth, $$ M_{\odot}=2\times 10^{30}\mathrm{Kg}, $$ the mass of the Sun and $$ r=1.5\times 10^{11}\mathrm{m}, $$ is the distance between Earth and the Sun. The latter defines what we call an astronomical unit AU. From Newton's second law we have then for the \( x \) direction $$ \frac{d^2x}{dt^2}=\frac{F_{x}}{M_E}, $$ and $$ \frac{d^2y}{dt^2}=\frac{F_{y}}{M_E}, $$ for the \( y \) direction.
Introducing \( x=r\cos{(\theta)} \), \( y=r\sin{(\theta)} \) and $$ r = \sqrt{x^2+y^2}, $$ we can rewrite $$ F_{x}=-\frac{GM_{\odot}M_E}{r^2}\cos{(\theta)}=-\frac{GM_{\odot}M_E}{r^3}x, $$ and $$ F_{y}=-\frac{GM_{\odot}M_E}{r^2}\sin{(\theta)}=-\frac{GM_{\odot}M_E}{r^3}y, $$ for the \( y \) direction.
We can rewrite these two equations $$ F_{x}=-\frac{GM_{\odot}M_E}{r^2}\cos{(\theta)}=-\frac{GM_{\odot}M_E}{r^3}x, $$ and $$ F_{y}=-\frac{GM_{\odot}M_E}{r^2}\sin{(\theta)}=-\frac{GM_{\odot}M_E}{r^3}y, $$ as four first-order coupled differential equations $$ \frac{dv_x}{dt}=-\frac{GM_{\odot}}{r^3}x, $$ $$ \frac{dx}{dt}=v_x, $$ $$ \frac{dv_y}{dt}=-\frac{GM_{\odot}}{r^3}y, $$ $$ \frac{dy}{dt}=v_y. $$
The four coupled differential equations $$ \frac{dv_x}{dt}=-\frac{GM_{\odot}}{r^3}x, $$ $$ \frac{dx}{dt}=v_x, $$ $$ \frac{dv_y}{dt}=-\frac{GM_{\odot}}{r^3}y, $$ $$ \frac{dy}{dt}=v_y, $$ can be turned into dimensionless equations (as we did in project 2) or we can introduce astronomical units with \( 1 \) AU = \( 1.5\times 10^{11} \).
Using the equations from circular motion (with \( r =1\mathrm{AU} \)) $$ \frac{M_E v^2}{r} = F = \frac{GM_{\odot}M_E}{r^2}, $$ we have $$ GM_{\odot}=v^2r, $$ and using that the velocity of Earth (assuming circular motion) is \( v = 2\pi r/\mathrm{yr}=2\pi\mathrm{AU}/\mathrm{yr} \), we have $$ GM_{\odot}= v^2r = 4\pi^2 \frac{(\mathrm{AU})^3}{\mathrm{yr}^2}. $$
The four coupled differential equations can then be discretized using Euler's method as (with step length \( h \)) $$ v_{x,i+1}=v_{x,i}-h\frac{4\pi^2}{r_i^3}x_i, $$ $$ x_{i+1}=x_i+hv_{x,i}, $$ $$ v_{y,i+1}=v_{y,i}-h\frac{4\pi^2}{r_i^3}y_i, $$ $$ y_{i+1}=y_i+hv_{y,i}, $$
It is rather straightforward to add a new planet, say Jupiter. Jupiter has mass $$ M_J=1.9\times 10^{27}\mathrm{kg}, $$ and distance to the Sun of \( 5.2 \) AU. The additional gravitational force the Earth feels from Jupiter in the \( x \)-direction is $$ F_{x}^{EJ}=-\frac{GM_JM_E}{r_{EJ}^3}(x_E-x_J), $$ where \( E \) stands for Earth, \( J \) for Jupiter, \( r_{EJ} \) is distance between Earth and Jupiter $$ r_{EJ} = \sqrt{(x_E-x_J)^2+(y_E-y_J)^2}, $$ and \( x_E \) and \( y_E \) are the \( x \) and \( y \) coordinates of Earth, respectively, and \( x_J \) and \( y_J \) are the \( x \) and \( y \) coordinates of Jupiter, respectively. The \( x \)-component of the velocity of Earth changes thus to $$ \frac{dv_x^E}{dt}=-\frac{GM_{\odot}}{r^3}x_E-\frac{GM_J}{r_{EJ}^3}(x_E-x_J). $$
We can rewrite $$ \frac{dv_x^E}{dt}=-\frac{GM_{\odot}}{r^3}x_E-\frac{GM_J}{r_{EJ}^3}(x_E-x_J). $$ to $$ \frac{dv_x^E}{dt}=-\frac{4\pi^2}{r^3}x_E-\frac{4\pi^2M_J/M_{\odot}}{r_{EJ}^3}(x_E-x_J), $$ where we used $$ GM_J = GM_{\odot}\left(\frac{M_J}{M_{\odot}}\right)=4\pi^2 \frac{M_J}{M_{\odot}}. $$ Similarly, for the velocity in \( y \)-direction we have $$ \frac{dv_y^E}{dt}=-\frac{4\pi^2}{r^3}y_E-\frac{4\pi^2M_J/M_{\odot}}{r_{EJ}^3}(y_E-y_J). $$ Similar expressions apply for Jupiter. The equations for \( x \) and \( y \) derivatives are unchanged. This equations are similar for all other planets and as we will see later, it will be convenient to object orient this part when we program the full solar system.
NASA has an excellent site at http://ssd.jpl.nasa.gov/horizons.cgi#top. From there you can extract initial conditions in order to start your differential equation solver. At the above website you need to change from OBSERVER to VECTOR and then write in the planet you are interested in. The generated data contain the \( x \), \( y \) and \( z \) values as well as their corresponding velocities. The velocities are in units of AU per day. Alternatively they can be obtained in terms of km and km/s.
For the first simple system involving the Earth and the Sun, you could just initialize the position with say \( x=1 \) AU and \( y=0 \) AU.
When building up a numerical project there are several elements you should think of, amongst these we take the liberty of mentioning the following.
The simplest possible step is to code the Earth-Sun system using Euler's method in two dimensions. The four coupled differential equations can then be discretized using Euler's method as (with step length \( h \)) are $$ v_{x,i+1}=v_{x,i}-h\frac{4\pi^2}{r_i^3}x_i, $$ $$ x_{i+1}=x_i+hv_{x,i}, $$ $$ v_{y,i+1}=v_{y,i}-h\frac{4\pi^2}{r_i^3}y_i, $$ $$ y_{i+1}=y_i+hv_{y,i}, $$
/*
Euler's method for the Eart-Sun system, simplest possible code
*/
#include <cmath>
#include <iostream>
#include <fstream>
#include <iomanip>
using namespace std;
// output file as global variable
ofstream ofile;
// function declarations
void output( double, double, double, double, double);
int main(int argc, char* argv[])
{
// declarations of variables
char *outfilename;
// Read in output file, abort if there are too few command-line arguments
if( argc <= 1 ){
cout << "Bad Usage: " << argv[0] <<
" read also output file on same line" << endl;
// exit(1);
}
else{
outfilename=argv[1];
}
ofile.open(outfilename);
int NumberofSteps = 1000;
double FinalTime = 1.0;
double Step = FinalTime/((double) NumberofSteps);
double time = 0.0;
// Initial values x = 1.0 AU and vy = 2*pi
double pi = acos(-1.0);
double FourPi2 = 4*pi*pi;
double x = 1.0; double y = 0.0; double vx = 0.0; double vy = 2.0*pi;
double r = sqrt(x*x+y*y);
// now we start solving the differential equations using Euler's method
while (time <= FinalTime){
x += Step*vx;
y += Step*vy;
vx -= Step*FourPi2*x/(r*r*r);
vy -= Step*FourPi2*y/(r*r*r);
r = sqrt(x*x+y*y);
time += Step;
output(time, x, y, vx, vy); // write to file
}
ofile.close(); // close output file
return 0;
} // End of main function
// function to write out the final results
void output(double time, double x, double y, double vx, double vy)
{
ofile << setiosflags(ios::showpoint | ios::uppercase);
ofile << setw(15) << setprecision(8) << time;
ofile << setw(15) << setprecision(8) << x;
ofile << setw(15) << setprecision(8) << y;
ofile << setw(15) << setprecision(8) << vx;
ofile << setw(15) << setprecision(8) << vy << endl;
} // end of function output
In case the function to integrate varies slowly or fast in different integration domains, adaptive methods are normally used. One strategy is always to decrease the step size. As we have seen earlier, this leads to more CPU cycles and may lead to loss or numerical precision. An alternative is to use higher-order RK methods for example. However, this leads again to more cycles, furthermore, there is no guarantee that higher-order leads to an improved error.
Assume the exact result is \( \tilde{x} \) and that we are using an RKM method. Suppose we run two calculations, one with \( h \) (called \( x_1 \)) and one with \( h/2 \) (called \( x_2 \)). Then $$ \tilde{x}=x_1+Ch^{M+1}+O(h^{M+2}), $$ and $$ \tilde{x}=x_2+2C(h/2)^{M+1}+O(h^{M+2}), $$ with \( C \) a constant. Note that we calculate two halves in the last equation. We get then $$ |x_1-x_2| = Ch^{M+1}(1-\frac{1}{2^M}). $$ yielding $$ C=\frac{|x_1-x_2|}{(1-2^{-M})h^{M+1}}. $$ We rewrite $$ \tilde{x}=x_2+\epsilon+O((h)^{M+2}), $$ with $$ \epsilon = \frac{|x_1-x_2|}{2^M-1}. $$
With RK4 the expressions become $$ \tilde{x}=x_2+\epsilon+O((h)^{6}), $$ with $$ \epsilon = \frac{|x_1-x_2|}{15}. $$ The estimate is one order higher than the original RK4. But this method is normally rather inefficient since it requires a lot of computations. We solve typically the equation three times at each time step. However, we can compare the estimate \( \epsilon \) with some by us given accuracy \( \xi \). We can then ask the question: what is, with a given \( x_j \) and \( t_j \), the largest possible step size \( \tilde{h} \) that leads to a truncation error below \( \xi \)? We want $$ C\tilde{h} \le \xi, $$ which leads to $$ \left(\frac{\tilde{h}}{h}\right)^{M+1}\frac{|x_1-x_2|}{(1-2^{-M})}\le \xi, $$ meaning that $$ \tilde{h}=h\left(\frac{\xi}{\epsilon}\right)^{1+1/M}. $$
With $$ \tilde{h}=h\left(\frac{\xi}{\epsilon}\right)^{1+1/M}. $$ we can design the following algorithm: