My initial intent was to write a series of posts detailing some of the functionality of the Sundials package released by Lawrence Livermore National Laboratory (https://computation.llnl.gov/projects/sundials). Sundials is a package for numerically solving ordinary differential equations and differential algebraic equations. It also comes with a library for the solution of nonlinear algebraic equations. We will go into detail on this at a later date.
I've used Sundials once or twice, but some current projects I'm working on require the speed and horsepower of C++ as opposed to MATLAB or Python. I was making more of effort to familiarize myself with the full capabilities of the package. To that end I was working through a couple of simple problems: a projectile motion problem with air resistance, and a damped driven harmonic oscillator. The driven oscillator problem is nonlinear and is often used as an introduction to chaos in courses on computational physics. Small changes in initial conditions or other parameters can lead to vastly different behaviors.
I started on the oscillator problem and grabbed an old book from my undergraduate days, Computational Physics by Nicholas J. Giordano. I noticed something when reworking the example presented by the author. As mentioned, the problem is a driven pendulum consisting of a mass hanging from a rigid rod of length $L$. The angle between the rod and the vertical is denoted by $\theta$. This is shown below in Fig. 1.
 |
| Figure 1: Schematic of pendulum |
The author considers the case where the pendulum is driven by a periodic force with amplitude $F_d$ and frequency $\Omega$,
$$F_\mbox{drive} = F_d \sin (\Omega_d t).$$
We will also include a damping force proportional to the angular velocity,
$$F_{\mbox{damping}} = -q \omega = -q \frac{d\theta}{dt},$$
where $q$ is the constant of proportionality and $\omega$ is the pendulum's angular velocity.
The equation of motion of the pendulum is given by,
\begin{equation}
\frac{d^2 \theta}{dt^2} = -\frac{g}{L}\sin{(\theta)} - q \frac{d\theta}{dt} + F_d\sin{(\Omega_D t)}
\label{Main_de}
\end{equation}
subject to the initial conditions for our starting angle and angular velocity, $\theta(0) = \theta_0$ and $d\theta(0)/dt = \omega_0$.
We can rewrite Eq. \ref{Main_de} as two coupled first-order differential equations,
$$\frac{d \omega}{dt} = -\frac{g}{L}\sin{(\theta)} - q \frac{d\theta}{dt} + F_d\sin{(\Omega_D t)} $$ and
$$\frac{d\theta}{dt} = \omega.$$
The author solves these equations using an Euler-Cromer algorithm where we step a small time interval $\Delta t$, and calculate approximate solutions at time $t_{i+1} = t_i + \Delta t$ as follows,
$$\omega_{i+1} = \omega_i - \frac{g}{L}\sin(\theta_i)\Delta t - q\omega_i \Delta t + F_d\sin(\Omega_d t) \Delta t$$
$$\theta_{i+1} = \theta_i + \omega_{i+1}\Delta t.$$
Below is C++ code (a translation of the code in the textbook) implementing the above algorithm for the following model parameters.
- g = 9.8 $\mbox{m}/s^2$
- L = 9.8 m
- $F_d = 1.2 \;\mbox{s}^{-2}$
- $q_d = 0.5 \;\mbox{s}^{-1}$
- $\Omega_d = 2/3 \;\mbox{s}^{-1}$
- $\theta(0) = 0.2$
- $\omega(0) = 0\;\mbox{s}^{-1}$
#include <iostream> #include <math.h> using namespace std; int main() { double g = 9.8; double L = 9.8; double F_d = 1.2; double q_d = 0.5; double Omega_d = 2.0/3.0; double period = 2 * M_PI * sqrt(g/L); int i = 0; double t = 0; double dt = 0.04; double theta = 0.2; double omega = 0.0; while (t <= 10 * period) { i = i + 1; t = t + dt; omega = omega - g/L * sin(theta) * dt - q_d * omega * dt + F_d * sin( Omega_d * t) * dt; theta = theta + omega * dt; cout << t << "\t" << "\t" << theta << "\t" << omega << endl; } return 0; }
Fig. 2 shows the results of this integration.
 |
| Figure 2. Results of an Euler-Cromer integration of the pendulum's equations of motion with a time step of 0.04 seconds. |
I mentioned this came about in an attempt to write an article detailing some of the functionality of Sundials. For the sake of argument, let's take the output of Sundials as correct and overlay those results on the data obtained from Euler-Cromer.
 |
| Figure 3. Euler-Cromer integration with the results from Sundials overlay ed. The two methods agree up to around 30 seconds, then begin to drastically diverge. |
OK, there is an obvious discrepancy here. After starting out the same, the two different simulations diverge considerably as time goes on. Is it the algorithm itself? Well, let's consider the case with no driving force and dissipation. Let's also take the initial displacement $\theta$ to be small such that we can make use of the approximation $\sin(\theta) \approx \theta$. In other words, we want to solve the equation,
\begin{equation}
\frac{d^2 \theta}{dt^2} = -\frac{g}{L}\theta.
\label{simple_de}
\end{equation}
This does have an analytical solution,
$$\theta(t) = \theta_0 \sin\left(\sqrt{\frac{g}{L}} t + \frac{\pi}{2}\right).$$
Let's rerun our code with the driving force and dissipation set to zero and plot Eq. \ref{simple_de} on top of it.
 |
| Figure 4: Comparison between analytical results fro a simple harmonic oscillator and the numerical solution using the Euler-Cromer algorithm. |
Well, OK. Those two curves are identical so the issue probably isn't with the code itself. The issue is the time step used in the Euler-Cromer code. As I mentioned earlier, the pendulum problem is often used as an introduction to chaotic behavior-- slight differences in parameters or initial conditions can lead to radically different behavior. We can see that by running the sundials code for slightly differnt initial starting angles. In the plot below, we start with our original initial position and solve for three cases where we increment the initial starting angle by 2%. This works out to be slightly less than 0.25 degrees each time. That isn't a big change, so we shouldn't see much of a change. Indeed, we don't see that big of a difference between the plots until roughly 35 seconds as the four curves begin to diverge considerably.
 |
| Figure 5: Solution to the damped driven harmonic oscillator problem for four slightly different starting angles. The solutions are similar until around 35 seconds at which point they begin to differ significantly. |
The Euler-Cromer algorithm isn't particularly accurate, so errors in the results at each time step can accumulate faster than in other, more accurate integration schemes. This isn't really important in linear problems-- the simple harmonic oscillator plotted above is just fine-- but when we work with the nonlinear problems those small errors change the predicted results entirely. We can confirm this is indeed what is happening by reducing the time step in the Euler-Cromer integration. In Figure 6, we take this time step to be 0.0001 seconds. As can be seen, this now matches well with the results out of Sundials, although some discrepancy is creeping in towards the end of the simulation.
 |
| Figure 6: Comparison between Sundials and Euler-Cromer algorithm with the time step in Euler-Cromer set to 0.0001 S |
 |
| Figure 7: Plot of solutions using relative tolerances of 1e-3 and 1e-6 as well as the solution as calculated by Sundials. For clarity the density of plotted points is lower than in previous plots. |
The bottom line: Even simple looking systems can produce unexpected results and be very sensitive to the algorithms (and internal settings of the algorithms) used to solve them. Don't just blindly trust the output. Always double check and build in various means of sanity checking your results.
No comments:
Post a Comment