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 realtime 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 readytouse 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.
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 righthand sides. The parametric equations must be in the Auxiliary variables box, and just one trivial ODE:
Then you can compile and graph this problem using all the advanced graphic features of this system.
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.
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.
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 PseudoSphere PseudoSphere.scr .
In cases when we need to outline a stripe, we can use a sine wave with a short period and the amplitude equal to halfwidth 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 "onesided" 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.
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 Doublependulum 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 Doublependulum (to derive them is a challenge and they are quite cumbersome).
Slightly digressing, make a note about the pendulumlike 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 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).
A large number of simulations here are dedicated to various cases of the Newtonian nbody 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 nonplanar 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 i_{15} used in the equations (this constant is neg1p5 in autogenerated nbody problems). Presuming that the central force is proportional to r^{k}, i_{15} = (k – 1)/2. For the Newtonian central force, k = –2 so that i_{15} = –1.5. For nonNewtonian central force i_{15} differs from –1.5. In particular, for a linear spring k = 1 so that i_{15} = 0.
The constant i_{15} appears as an exponent in the equations. Therefore, if i_{15} = 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 2body 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 r^{k}. Explain students what is the central force generated by a pointed source: the force proportional to r^{k}, 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 i_{15} , demonstrate, that elliptic trajectories are possible only for the special values k = –2 and k = 1. For all other values of k ellipselike 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 nbody 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 3body Lagrange case either from the Demo menu, or load scripts 3Bodies2D.scr, 20BodiesCirc.scr, 20BodiesEll.scr, or generate any desired plane nbody 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 threebody 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 5degree 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 5degree polynomials associated with the Euler case.
7. The Lagrange points. Here is a remarkable problem for obtaining of the socalled 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 guessvalue 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 ellipselike 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 toruslike shape – though not exactly a torus: 4BodiesTorus.scr .
In recent decades scientist have discovered never expected regular and periodic solutions with quite remarkable features in the threebody problem. Among the classical solution, the only known periodic solutions of the threebody 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 nonelliptic 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 subclass of nonelliptic periodic orbits was discovered by Xiaoming LI and Shijun LIAO: the Freefall 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 nonplanar choreographies for more than three bodies.
Alexander Gofen
April 2024