Ideas for courses in applied ODEs

Based on the preloaded simulations


This software comes with a wide list of preloaded simulations – i.e. script files running the integration and graphing of particular problems, and then allowing to Play the solutions as real-time animations. Such animations or simulations for problems taught in ODE courses increase motivation and interest in students. Moreover, the software allows to change the parameters and experiment with the studied processes as if in a virtual lab, or Exploratorium.  

Here is a special page Exploratorium in the Taylor Center offering several ready-to-use textbooks supported by the respective simulations.  However, the list of preloaded simulations is much wider containing over 200 simulation not supported with any didactical material. This is merely an annotated list of the script files in various topics of applied ODEs, so that it is a matter of future to develop the respective textbooks for those simulations analogous to those in the Exploratorium.

With that in mind, I compiled a list of possible syllabi for various courses grouped by topics, to be developed into textbooks for teaching the theory of those 200 preloaded simulations.  

When you load script files of simulations or open ODE files coming with this software, all the respective equations may be viewed and studied in the respective boxes in the front panel under the Tab Equations setting.



Analytic geometry


As a side application, this software may be used merely for drawing 2D or 3D curves defined by their parametric equations  

x = f(t),

y = g(t)

z = h(t)


where functions f , g, and h are expressions allowed in the right-hand sides. The parametric equations must be in the Auxiliary variables box, and just one trivial ODE:



x = f(t)

y = g(t)

z = h(t)

t = t0       

t' = 1



Then you can compile and graph this problem using all the advanced graphic features of this system.

Drawing 2D and 3D curves

Using the parametric equations, you can visualize the curves that you study (or suggest your students to do it): say the Corny spiral, or Curly spirals (CornuSpiral.scr, CurlySpiral.scr) both available in the Demo menu. See the notes of theory for those curves here.


Drawing 3D curves which outline surfaces

Using the parametric equations for curves, you can also visualize 3D surfaces outlining them either with one curve, or with a family of curves, as it will be explained below. Outlining surfaces with 3D curves is a good method of visualization and exercise in connection with the equations of surfaces, and this method also allows to "see through" the surfaces.

Using a helix as an outline

1.      Visualize two linked tori as in the Demo menu (Tori.scr).

2.      Visualize the Klein bottle (ScrewLineKleinBottle.scr) whose equations are here (scroll down).

3.      Visualize the Pseudo-Sphere PseudoSphere.scr .


Using a sine wave as an outline

In cases when we need to outline a stripe, we can use a sine wave with a short period and the amplitude equal to half-width of the stripe. For example, so we can visualize the Möbius surface (MobiusSineOutline.scr).

In all samples visualizing Möbius surface you can modify its number of twists changing the constant n (for the simplest "one-sided" Möbius surface n = 0.5).

The method of outlining by helix and sine wave is convenient because it requires to specify only one curve to outline the entire surface, while families of curves (below) require to specify sets of curves resorting to the feature Arrays of IVPs in this software.


 Using families of curves as outline

Traditionally, surfaces are outlined using two families of curves corresponding to the curve coordinates (u, v) on the surface:

x = f(u, v)

y = g(u, v)

z = h(u, v)

Here are a few such samples.

When the desired number of outlining curves is small, their equations may be entered manually. Such are Möbius samples Mobius.scr and MobiusLarge.scr .

However, this software also provides mechanism for cloning one IVP into an array of IVPs with different parameters, and this approach was applied in samples ToriArray.ode, MobiusArray.ode, and   KleinBottleForArray.ode . These ODE files include special lines for generating arrays. Before compiling them, first you have to go to …

·        Create array of IVPs

·        Clone excluding independent variable

·        Click "Unfold" button.


Then Compile, click 3D Check box, click the button {x1, y1} >>, and the button Graph, which will display the 3D anaglyphic stereo picture of the family of the curves.

You can skip these steps loading the full scripts  FullTori.scr, FullKleinBottle.scr .



First, illustrate the classical mathematical pendulum loading Pendulum2D.scr . In it the constant Te0 stands for the initial angle of declination from the vertical, and Vte0 stands for the initial angular velocity. In the preloaded sample Te0 < π and the initial velocity is zero – as it is expected from a pendulum. You can set  Te0 = π and watch how long it will take the rounding errors of the computation to push the pendulum off this point of unstable equilibrium (set the Play time about 15 s). Then experiment with a nonzero velocity demonstrating that under certain conditions the pendulum spins around the pivot rather than swings.

Then consider a couple of examples of 3D pendulum  PendulumApple.scr   and  PendulumFlower.scr . Here you can play and explain the meaning of two initial angular velocities also, demonstrating the conditions when the pendulum spins rather than swings.  

Then derive the ODEs and demonstrate a plane Double-pendulum loading DoubePendulum.scr . In the preloaded sample the parameters are such, that the pendulum demonstrates chaotic motion. Suggest your students to figure out the initial values which generate a few types of regular motion, and to try them in the simulation. Finally, suggest the students to make a search if they can find the ODEs for 3D Double-pendulum (to derive them is a challenge and they are quite cumbersome).

Slightly digressing, make a note about the pendulum-like behavior of the Heavy top (considered in the three cases of a rigid body motion in the Exploratorium). Though normally it is presumed that a Heavy top rests on a pivotal point above the horizontal surface, its ODEs also model the case as if it hung at the pivotal point below the horizontal surface swing about this point as a pendulum. There is the analysis and the supporting simulations in the textbook demonstrating how the big amount of spin of the heavy top prevents its falling, while a small amount or no spin makes it behave as a pendulum.

The Lorenz attractor


The Demo menu offers a demonstration of the widely known Lorenz attractor, for which teachers may easily find the supporting theory, and then play it varying the parameters. 


An amazing process of Autocatalators is explained in the paper by the late Borrelli and Coleman "Pitfalls and Pluses in Using Numerical Software to Solve Differential Equations", example 3, page 11 (2009).  See the simulation illustrating this reaction loading script WineGlassx(t)y(t).scr for dependency on t. The "nickname" Wine Glass becomes obvious when you view the simulation WineGlass.scr as y(x).  

The Newtonian n-body problem


A large number of simulations here are dedicated to various cases of the Newtonian n-body problem containing their own set of the ODEs for every particular case.

Besides them, this software provides a menu item for automatic generation of a set of Newtonian ODEs for any number of bodies n [2..99] for planar or non-planar cases. This automation sets the initial values as for the Lagrange setting (elliptic or circular), but you then can change these values according with your needs.  


In all problems for Newtonian motion, there is a special constant  i15  used in the equations

(this constant is neg1p5 in auto-generated n-body problems).

Presuming that the central force is proportional to  rk,  i15 = (k – 1)/2.

For the Newtonian central force,  k = –2  so that  i15 = –1.5.

For non-Newtonian central force  i15  differs from –1.5.

In particular, for a linear spring   k = 1 so that  i15 = 0.


The constant i15 appears as an exponent in the equations. Therefore, if  i15 = 0, the respective factors turn into 1 becoming redundant  so that the system may be simplified manually by removing those factors. However, it will compile and run properly also without simplifications.



1.           Infinitely small masses. Explain students a remarkable feature of the Newtonian equations in cases when one or several masses are infinitely small: for example, a mass of an artificial satellite vs. the mass of Earth. When any of the masses approaches zero, the Newtonian equations admit " its limit format" allowing to replace any occurrence of such mass with zero simplifying the system. Therefore, here are many simulations, where some masses in the Constants are set to zeros, or the equations were manually simplified in such cases. This must not confuse the students, that in the Newtonian mathematical model even zero masses are attracted by gravity. Such are simulations 2BodiesProbe.scr, 2plusProbe.ode exemplifying the simplest case of 2-body problem. In them you can demonstrate how various values of the initial velocity of a probe lead to orbits representing all types of conic sections. Draw students' attention also to accelerations and deceleration of the probe during motion along the orbit in accordance with the Kepler's laws.

2.           The central force proportional to  rk. Explain students what is the central force generated by a pointed source: the force proportional to  rk, and that the two fundamental forces of nature (gravitation and electric attraction) correspond to k = –2. Explain what is the conservative field. Then prove the remarkable fact, that the central force makes conservative field only for the values k = –2 (gravity or electrostatic) and  k = 1 (a linear spring – simulation NonNewton0.scr). The trajectories for both kinds are conic sections – demonstrate it with the simulations in item 1. Draw attention to the fact, that a mass exerting gravity is in a focus of an elliptic orbit, while for linear spring its pivot point is in the center of ellipse. Then, varying the value  i15 , demonstrate, that elliptic trajectories are possible only for the special values k = –2 and  k = 1. For all other values of k ellipse-like shapes drift (making a precession). Explain that Mercury's orbit does make such a precession, though not because k ≠ –2, but because the gravity of the Sun is so strong near Mercury, that the Newtonian mechanics must be replaced with the more accurate Einstein's General Relativity theory.

3.           The Lagrange case. Explain, that the n-body problem may be solved in elementary functions only for the number of masses n = 2. Research for n > 2 has been conducted for about 300 years since Newton. The solutions in closed format for cases of n > 2 were found only for a few special settings considered in the following items here. The simplest (and most intuitive of them) is the Lagrange case for n equal masses placed at the vertices of a regular polygon – see the definition and the didactic text here. Play the 3-body Lagrange case either from the Demo menu, or load scripts 3Bodies2D.scr, 20BodiesCirc.scr, 20BodiesEll.scr, or generate any desired plane n-body problem as provided in the menu. Draw attention to the fact, that the Lagrange orbits are unstable: the more unstable the more elongated the ellipses are. Use the available simulations of the Lagrange setting to demonstrate how the orbits become chaotic after several turns. The demonstrate the slightly disturbed Lagrange cases in 2D and 3D, which promptly become chaotic. Explain, that some of chaotic orbits, after remaining bounded for certain time, end up with an engaged pair of bodies orbiting around each other and a single body escaping to infinity in opposite directions (like in the 2D case). However, it is not known what is the criterion for the initial conditions that a chaotic orbits stay bounded forever.    

4.           A provocative question. Pose for students the following provocative question. They have just seen the classical Lagrange case for plane regular polygons. Is it possible to generalize the plane Lagrange case into a spatial one placing equal masses at vertices of any of the five Platonic polyhedrons? Which initial velocities is it necessary to set in order that the bodies continue to move preserving the initial formation of the same Platonic polyhedron? For example, there is a script file 4BodiesTetra.scr setting some initial velocities. However, it is clearly seen that they do not preserve the formation. Well, it is proved in the didactic text here, that the spatial Lagrange case is impossible – except the trivial radial motion towards the collision – like in this simulation  4BodiesTetraCollision.scr . This impossibility follows from simple geometric reasoning given here.

5.           The Lagrange case is possible also for bodies of unequal masses, though only for n = 3. During the motion, the three bodies preserve the formation of equilateral triangle moving along ellipses of different sizes (or along the conic sections). Given the particular masses, the initial positions and velocities of the bodies must be computed in accordance with the special formulas given in the literature and used in the section Constants in two available simulations: elliptic and circular, accessible via the Demo menu. In the case of circular motion, the three bodies move along the concentric circles (preserving an equilateral triangle).

6.           The Euler case. The other classical case of known three-body solutions is that of Euler, when the three body of equal or unequal masses preserve the initial formation in line – both elliptic and circular simulations are accessible via the Demo menu. When it is the three masses which are given, it is necessary to solve a polynomial equation of degree 5 in order to obtain the appropriate initial positions and velocities. Otherwise, it is possible first to choose the initial positions, and then to obtain the appropriate masses avoiding solving the 5-degree polynomial – the way it is done in both simulations. There are also a few simulations in the subfolder EulerWorkshop dedicated to studying and graphing those special 5-degree polynomials associated with the Euler case.

7.           The Lagrange points. Here is a remarkable problem for obtaining of the so-called five Lagrange points having a special significance in space engineering explained below. According to the Kepler laws, planets, or other bodies orbiting the Sun, cannot stay on the same radius except for a short while, because their angular velocities diminish with their growing distance from the Sun. However, in the problem of 7 bodies where a pair of bodies have finite masses orbiting each other (say the Sun and Earth), while the 5 remaining bodies are massless, it is possible that 3 of those 5 massless bodies stay on the same diameter forever (though their positions are unstable so that orbital corrections will be needed from time to time). Here are two simulations illustrating the Lagrange points for two versions of the pairs of masses – both accessible via Demo menu. According to the theory, the initial values for the Lagrange points should be obtained via solution of a polynomial of degree 5. The subfolder LagrangePointsWorkshop contains three IVPs used for obtaining the solution of the degree 5 polynomial without iterations resorting to the special feature of this software. This feature allows to switch independent variables and integrate ODEs for the inverse function. In this case, when we make the independent variable the value of the polynomial, we integrate from some guess-value to zero. In both cases the simulations display circular motion of the 7 bodies where 5 of them stay on the same diameter. This solution for the theoretical Lagrange points is an idealization of the more complex real space engineering problem (say the Webb telescope), where the Earth moves along an ellipse-like orbit having the Moon as its satellite, and other planets affect the solution too.   

8.           The isosceles settings. Consider in comparison two special isosceles settings of three bodies, one of them oscillating: spatial setting called the Sitnikov case vs. the plane case when the bodies collide – script file PlaneCollision.scr . In it, the three bodies in the isosceles setting with zero initial velocities freely fall towards each other while the middle body oscillates along a straight line – until the two bodies collied. On the contrary, in the spatial Sitnikov arrangement, initially, a pair of equal masses in a plane are set into an elliptic motion, while the third body moves along a perpendicular erected from the center of the mass of the pair – simulations Slingshot.scr, SlingshotEscape.scr . and SlingshotZeroMass.scr . In the latter case with zero central mass the pair stays in a plane performing the elliptic movement, while the middle body oscillates vertically in an aperiodic manner. However, when the middle body has nonzero mass, it does affect the motion of the pair so that it is not more planar nor elliptic. The pair jumps up and down in accordance with aperiodic vertical oscillation of the middle body. Unlike the case of zero mass of the middle body not affecting the pair so that the smallest distance between the bodies of the pair always remains the same, the nonzero mass of the middle bodies attracts the pair and may make them to approach each other closer and closer – until the next passage of the central body accelerates it so much that it escapes to infinity while the pair also escapes in the opposite direction – see SlingshotEscape.scr .   

9.           Consider a setting of two pairs where the bodies in each pair orbit around each other, while the centers of masses of the pairs also orbit each other in a larger circle so that the orbits of these four bodies outline a torus-like shape – though not exactly a torus: 4BodiesTorus.scr .   


Recent discoveries


In recent decades scientist have discovered never expected regular and periodic solutions with quite remarkable features in the three-body problem. Among the classical solution, the only known periodic solutions of the three-body problem were the elliptic solutions of Lagrange and Euler, and it was obvious that randomly chosen initial values led to chaotic orbits. So more surprising were the discoveries in the recent decades of various classes of periodic orbits other than ellipses. For some of those discoveries listed below the simulations are included into the set of preloaded samples coming with this software.

1.           A wide class of non-elliptic periodic orbits was discovered by Ana Hudomal, and their simulations are available under Demo menu 3 bodies/Periodic (closed). There is also a demo for her more sophisticated relatively periodic orbit.

2.           Another wide sub-class of non-elliptic periodic orbits was discovered by Xiaoming LI and Shijun LIAO: the Free-fall periodic orbits. Their simulations are available under Demo menu 3 bodies/Periodic (Free fall): 30 cases of periodic free fall of equal masses. They are described also in the paper Dropping bodies – by Richard Montgomery, and in Exploratorium. Explain your students, that orbits of dropped 3 bodies are obvious only for equilateral setting (ending with triple collision), and for isosceles setting (ending with double collision). Therefore, the discovered initial settings of three equal masses are all irregular triangles. At that, for a random setting the orbits of free fall are typically chaotic.

3.           Choreography is a special periodic motion of three bodies one after the other along the same orbit – see the source here. A trivial case of Choreography of n equal masses is the circular Lagrange orbit of those masses. Nobody expected anything else until the discovery the "shape 8" choreography by Montgomery and Simo. The Demo provides an access to over 300 samples of different plane choreographies of three bodies. Now scientists discovered also non-planar choreographies for more than three bodies.  

Alexander Gofen 

April 2024