9.5 Simulation of Ordinary Differential Equations

9.5.1 Comparison of Euler, Heun and Runge-Kutta Methods

The interactive picture in Fig.9.7 leads to a simulation, which shows the three methods in parallel for the example of the exponential function. The initial value for x = 0 is y0 = 1 but can also chosen differently. The number intervals in the variable region can be chosen between 1 and 24.

There are 4 straight lines in the picture , that can be drawn and turned with the mouse. These allow to easily visualize the construction of the approximations.

For the initially shown rough resolution one clearly recognizes the different convergence quality of the methods and the large superiority of Runge-Kutta – with the eye its error can not be noticed any more.

The description of the simulation contains further details and suggestions for experiments. It also contains a description of the complete codes, which are in each case a few lines that are repeated in a loop once for each point of computation. The calculation happens so quickly, that one does not notice the time development in this example. For the commercial programs one can specify how many points must be calculated per minute in order to create in the resulting graphs the impression of a temporal sequence.

In practice nowadays one does not have to make the effort to write computation algorithms for the solution of ordinary differential equations, since they can simply be called in all numerical programs via specifying a name. However one should have understood once how this “witchcraft” actually comes to be.

Figure 9.7: Comparison of the convergence of Euler (blue),Heun(green) and Runge-Kutta(red) methods for y = y (exponential function, blue line). The green lines can be pulled with the mouse to make the construction of the methods for one interval by hand possible. With a slider the number of intervals in the constant variable region can be changed (number of points n = number of intervals + 1). With the second slider y0 is adjusted (in the picture y0 = 1).

In Fig.9.8 the relative error of the three methods discussed in Fig.9.7 is shown, i.e. for example (Euler - y0ex)y0ex. The ordinate region has been spread, such that the differences are more visible. For the small number of two to three points (one to two intervals) in the variable region even the small error of the Runge-Kutta method becomes noticeable.

In order to also rate the error for a larger resolution, a second window of the simulation shows the relative deviation at the end of the last interval with high accuracy.

Figure 9.8: Comparison of the relative error for four intervals in the variable region (blue Euler, green Heun, red Runge-Kutta). The points show the relative deviations from the analytical value of the exponential function, that is shown in blue. For the Runge-Kutta method the error at the end of the interval is shown in a number field on the right. The scale on the y-axis depends on the accuracy achieved via the Euler method.

9.5.2 Differential Equations of First Order

We are using a Runge-Kutta procedure, that is integrated into EJS to visualize explicit differential equations of first order. Implicit equations play a minor role in elementary physics. Their numerical solution ca be achieved via iterations, that are built into the computational algorithms.

In the graphs we use for the variable the symbol x and for the ordinate the symbol y.

The following interactive Figure 9.9 shows the graph of a transient process, that is defined via the differential equation, which is shown in the text field y. In this presentation the individual computation points are shown; via the option boxes one can switch to a line presentation.

The differential equation shown in a text field can be edited or entered from scratch, such that you can investigate arbitrary explicit differential equations of first order. The speed of the animation and the accuracy of the calculation that is related to it can be varied with the slider for the step width.

Figure 9.9: Animated solution of differential equations of first order. The picture shows a convergent transient process. The range of the variable x, the initial value y0 and the step width for the calculation can be chosen. One can select either a point presentation or a smoothing line presentation. With the option boxes one can select a new calculation or a superposition of calculations with different settings.

The option box allows for the selection among a number of elementary differential equations, that are preset with initial value, for example

In the last two cases the solution of the differential equation is reduced to the normal integration process, since it does not contain y, and therefore these differential equations have as solutions the anti-derivatives of C and Cx.

The examples are classified according to the following characteristics:

Figure 9.10: Phase-space y versus y of differential equations of first order for the example of Fig.9.9. The green point shows during the animation, that starts at (1,1), the current computation point.

In Fig.9.10 a phase space projection y versus y is shown. This shows the character of the differential equation and of its solution , here oscillating convergent, quite clearly. The green point designates the current endpoint of the calculation. In this examples y converges against a finite value, while y converges to zero. en

The initial value y0 and the initial abscissa x0 can be chosen at will. The formulas are editable, such that you may enter arbitrary analytic functions and study them.

Multiple runs can be organized with the switches in order to compare the curves for different initial value, initial abscissae or differential equations. The passive picture in Fig.9.11 shows a simple example for constant acceleration a with the differential equation y = at. Here the acceleration a is increased over eleven steps from 0 to 10. The initial value stays at y0 = 5.

Figure 9.11: Family of solutions from the simulation of Fig.9.9 for the preset example y = ax (constant acceleration a) with parameter a = 0,1,,10 and y0 = 5. Interpretation: x= time(t) ,y distance, y0 start position.

The description pages of the simulation contain further details and many suggestions for experiments.

During the initial work with the simulation one is often surprised by the totally unexpected results when entering a specific equation for y or even only changing the value of a parameter in the equation. One is used to have a mental picture of dependencies of the form y(x), but this is not the case for y = f(x,y), if one is not familiar with that.

Thorough experimenting with this simulation and the following examples for the differential equations of second order is therefore necessary to obtain a thorough understanding of the relationships that are described by the differential equations.

These examples visualize, that a single differential equation defines a relationship, which contains an unlimited number of specific solutions. The initial values fix a particular solution from the family of possible solutions. The parameter that determines the specific solution does not have to be the initial value of the solution. One can also demand, that the value of the function must be y1 at a later time t1. For the numerical solution one then solves the differential equation starting from y1 firstly in the direction of increasing t > t1 and then in the direction of decreasing t < t1.

For a differential equation of first order the family of solutions has one parameter, for differential equations of second order the family of solutions has two parameters (see following subsection).

9.5.3 Differential Equations of second order

Numerous relationships in physics are described by differential equations of second order. In addition to the acceleration ( 2. derivative) they allow also to take velocity dependent interactions into account ( 1. derivative) that include friction processes. They also cover all undamped, purely periodic functions as special cases. Including damping makes realistic models of pendulums and oscillators possible.

The elementary functions described by differential equations of first order are covered by analogous differential equations of second order. We discussed already the differential equation of a similar structure will then contain additional functions. . The solutions of differential equations of second order constitute a 2 parameter family.

Only with two initial values y0 and y0 for the start value x0 of the variable a specific solution is fixed. Thus a single initial value still allows a whole 1-parameter family of solutions.

Among the explicit differential equations a very important one is the very simple one whose solutions are the trigonometric functions.

y = -y or y + y = 0 y = sint y = cost y = -sint = -y y = cost y = -sint y = -cost = -y y = eit = cost + isint y = ieit y = i2eit = -eit = -y

It describes many oscillation processes.

As already discussed at the end of section 2 of this chapter, one reduces this differential equation for the numerical solution to a system of two coupled differential equations of first order.

in general y = f(y,y,x) 1. definition: y(x) y1(x) 2. definition: y = y1 y2 y2 = y1 = y = f(y1,y2,x)

The first equation defines the first of the new functions by the original function. the second equation defines the second new function as derivative of the first one. The original differential equation connects y2 with y1,y2 and x . From the solution for y1,y2 one recovers y and y.

Thus the two coupled differential equations of first order  (a) y1 = y2 (b)y2 = f(y1,y2,x)for the two functions y1,y2  are equivalent to the single differential equations of second order of for y(x) y = f(y,y,x)

special case:  y = -y becomes via y y1and y = y1 = y2 the system of two differential equations for y1 and  y2  (a)y1 = y2; (b) y2 = -y1

The steps by which the two equations are solved for subsequent points have to be nested in a suitable way. Using equation (a) one first calculates an approximation for the derivative, which is then substituted in equation (b) instead of the formally required derivative. In practice this algorithm is contained in all popular numerical programs. We again use for our examples a EJS-simulation, for which we only add equation (a) in an additional line. (As designation for the first derivative we use in the formula field “yprime” since Java cannot understand “y”.)

For a differential equation of higher order this method would have to be repeated for every further order and chained in an equivalent manner. Differential equations of higher order do however play no major role in physics.

The following interactive picture in Fig.9.12 leads to a simulation for differential equations of second order. It shows an exponentially damped periodical oscillation. In the differential equation y = -y- 0.2y shown in the text box the first term - y is responsible for generating a periodic function and the second term - y for the exponential decay, as familiar from the differential equations of first order. The factor 0.2 determines the speed of decay.

The control of this simulation is quite similar to the case of differential equations of first order; only control elements for the second initial value y are added. In the selection box differential equations and initial values for the following functions are preset:

All parameters can be changed. In the text-field the differential equation can be changed or a totally new one can be entered, such that you may investigate arbitrary differential equations of second order using this simulation.

Figure 9.12: Animated simulation of the solution of differential equations of second order. Example: damped oscillation with the initial values y0 = 0 and y0 = 1. The arrow shows the two initial values for value and derivative. Variable range, initial values, step width and presentation type can be set. In addition the phase space diagrams can be shown in a 2D or a 3D presentation.

With two switches a 2D-presentation y(y) and y(y) and/or a rotating 3D-presentation y(y,y) of the phase spaces can be chosen. This window of the simulation is shown in Fig.9.13

Figure 9.13: Animated phase spaces for the differential equations of Fig.9.12. The left windows plots y(y) in blue and y(y) in red. The respective end points are highlighted and marked in colour. In the right windows y, y and y are mapped to the three space axes. The red curve thus represents the total differential equation y = f(y,y). The 3D projection is rotated in the animation.

The two-dimensional phase diagram now shows two curves y(y) in red and y(y) in blue. In this example one recognizes the damped transient process as double exponential spiral.

The three-dimensional phase diagram shows y = f(y,y) as a plane spiral in phase space. Its rotation during the animation increases the spatial impression.

The description pages of the simulation file contain further details and suggestions for experiments.

9.5.4 Differential Equations for Oscillators and Gravity Pendulum

The differential equations of second order discussed in section 9.5.3 describes among other system all possible kinds of oscillators, including also the classical mathematical gravity pendulum (called mathematical, because it treats the pendulum as a mass point on a mass-less stiff rod in abstraction from its real construction). For these cases the differential equations and initial conditions of the following simulation are pre-formulated, which otherwise is very similar to the previous one.

The interactive picture in Fig.9.14 shows the example of a damped oscillator, which initially oscillates in its eigen frequency until x = 30, when an external force at double the frequency is added to the system. One sees the transition from the free oscillation to the forced oscillation at double the frequency including interferences. The free oscillation finally decays away totally. The driven oscillation remains with double the frequency and a constant amplitude.

The corresponding phase space curve in Fig.9.15 is quite confusing as a static picture. If one however observes the dynamic flow, one recognizes the different transitions quite easily.

Cleared of factors, that only scale the graphics and are needed for the formula to be recognized (yprime instead of y) the differential equation reads: y = -y - y + sin2xstep(x - 30).

Figure 9.14: In the picture the solution of an oscillation equation with damping, which is driven by an external force at double the eigen frequency and supplied with energy. The differential equation and all parameters can be changed.

The term - y produces a periodic oscillation with period 2π, the term - y an exponential damping and the term sin2x a driving force with constant amplitude and the period π. The very useful step function switches at the given point in time x = 30 from 0 to 1. The damped oscillation of the free pendulum simply continues, while the periodic driving force is added at this time.

In the phase space picture shown in Fig.9.15 one also recognizes the transition between the two kinds of oscillations, from the initially free and damped oscillation (initial plane spiral) to the forced oscillation. After a sufficiently long time the free oscillation has been damped away and the oscillator moves periodically with constant amplitude a the frequency of the driving force.

Figure 9.15: Phase space plots for the oscillation equation of Fig.9.14. On the left projections y versus y in blue and y versus y in red, on the right y versus y and y. The picture shows the state shortly after adding the external driving force, on the left as lines, on the right as sequence of calculated points

The simulation contains the following pre-defined oscillators:

In addition for the gravitation pendulum as second pendulum ( full period of 2 s), the following situations are preset:

The plots of phase space curves for the gravity pendulum in the passive figure 9.16 show in the left window the situation for a deflection of 5.7 degree, for which the oscillation is still practically sinusoidal (red curve y -y) and in the right window for a deflection of 57 degree, for which the oscillation deviates quite clearly from it. Thus the blue curve is therefore for the small deflection a circle which is traversed with constant angular velocity (please note the different scales on the axes!). For the large deflection on recognizes in the animated simulation the extended time spent in the vicinity of the turning points.

The formulas and initial values can again be changed. In the vicinity of the unstable equilibrium (deflection of π) the solutions become extremely sensitive to the initial values, but also to the accuracy of the computation, that can be adjusted with the step width slider.

The description pages contain again exact details and hints for experiments.

Figure 9.16: Phase space curves of y versus y in blue and of y versus y in red for the pendulum example of the simulation in Fig.9.14: shown on the left for small and on the right for large deflections. Please note the different scales on the axes, especially the ordinate scale. The red line in the left window shows the negative linear relationship between acceleration and angle of deflection. In the right window one recognizes the large non-linearity for large deflections. Therefore only pendulums with deflections of a few degrees can be used for accurate clocks

9.5.5 Character of Ordinary Linear Differential Equations

From the experimental analysis of different explicit linear differential equations of second order we can draw a few general conclusions:

the following term in the differential equation means respectively: y = -y periodic function with period 2π y = -a2y periodic function with period 2πa y = y exponential growth with x y = yexponential change with x

y = const. constant acceleration  y = 0 constant velocity (0 acceleration) y = f(x) x-dependent driving force, characterized by f(x)  y = -yf(x) periodic oscillation , moderated by f(x)  y = -yg(x) exponential decay , moderated by g(x)

The points to which convergent or divergent solutions move in the phase diagram are referred to as point attractors and the closed target curves of periodic solutions are called periodic attractors.

9.5.6 Chaotic Solutions of Coupled Differential Equations


A new phenomenon appears, if three or more differential equations of first order are coupled and contain terms that are nonlinear in the variables. For certain parameter regions or regions of the initial values or even for all initial values their solutions show chaotic behaviour. This is especially attractive for oscillating systems that are characterized by differential equations of second order with the fundamental dependence y = -y....

Driven Double Pendulum

As first example we want to investigate in Fig.9.17 the simulation of a double pendulum, for which a second mathematical pendulum is fixed to the end of a primary mathematical pendulum (mathematical means here: the total mass is concentrated at the end of a stiff mass- and weightless pendulum rod).

The primary pendulum can be driven by a periodic force. The secondary pendulum is driven by the primary pendulum. Both are subject to gravity.

Each pendulum is described by an ordinary differential equation of second order, which corresponds to four differential equations of first order, and the differential equations are coupled, thus they also contain variables of the the corresponding other pendulum. It is now essential, that these differential equations are coupled by trigonometric functions and quadratic terms:

y1 = f1(y1,siny2,y2,sin(y2 - y1),y12,y22) y2 = f2(y1,siny1,y2,sin(y2 - y1),y12,y22)

The exact formulas are discussed in the description pages of the simulation.

The ratio of the pendulum lengths and the pendulum masses can be adjusted as well as the speed of the animation.

Figure 9.17: Chaotic movement of a driven double pendulum with adjustable length and masses (red curve). In the left picture the double pendulum is shown ( pivot point in green, mass point of the primary pendulum in blue, mass point of the secondary pendulum in yellow, vector of the external driving force as a blue arrow). In the right window the phase space projection angular velocity dϕdt versus angle of deflection ϕ is shown.

The indirect external driving force modulates the angular velocity of the primary pendulum with a sine function of adjustable frequency. The blue arrow shows direction and absolute value of the external driving force.

The red curve shows the orbit of the secondary pendulum, that means of the mass point at the end of the mass-less pendulum rod of length l2. It is possible to superimpose orbits for different initial values and thus to study at the same time the influence of small changes in the initial conditions on the long time behaviour.

In the right coordinate system of Fig.9.17 a plane phase space projection for the orbit of the primary pendulum is plotted. In addition a rotating presentation of the three-dimensional phase space y versus y versus y can be switched on (see Fig.9.18).

There obviously is no periodic attractor. One refers to a strange attractor, if the phase space orbits of the process described are limited to a certain region of the phase space, but do not become periodic, but show a fractal character and therefore cannot be described in an analytic closed form.

Together with adjusting the ratios of pendulum length and mass one obtains a rich spectrum of oscillations processes that happen chaotic but strictly deterministically.

Figure 9.18: 3D phase space diagrams for the double pendulum, on the left red for the end of the primary pendulum and on the right for the end of the double pendulum

The Reset button resets the simulation exactly (within the accuracy of the PC) to the same initial conditions. You may convince yourself, that the time development, which looks so confused is indeed repeated, and thus happens deterministically and not controlled by chance (this observation is achieved easily by calling the simulation twice and letting it run twice).

You may however also adjust the initial conditions for the position manually via pulling the yellow point. Then you will not achieve an exact reproduction, and two simulations running in parallel will soon run apart from each other . Thus the chaotic-deterministic characters is connected to an extreme dependence on the initial conditions .

The simulation description contains many suggestions for experiments.

For the double pendulum with its many non linear connections it will be barely possible to find a setting, that leads to a periodic solution. In other cases there are regions of chaotic next to regions of periodic behaviour.

Reflection of a Ball between sloping walls

For the second example shown in Fig.9.19 it is obvious, that there must also be periodic solutions.

Figure 9.19: Reflection of a ball between two sloping walls. In addition to chaotic orbits there also are periodic solutions, for which the ball jumps regularly back and forth on periodical orbits between the walls.

In this Ball in Wedge simulation a ball is reflected back and forth between two infinitely extended planes. For an initial orbit that starts symmetrically to the axis of symmetry and orthogonal to one of the surfaces the orbit is already closed after hitting both surfaces once. It can be suspected, that there are further periodic orbits with many reflections. In general however, the orbits are chaotic. The pitch of the surfaces, the position and the initial velocity of the ball can be adjusted via pulling with the mouse. The non-linearity of the connections lies here in the trigonometric functions used in the three coupled differential equations of first order.

The example also demonstrates the use of the Poincar section for the visualization of chaotic or periodic orbits and its use for the determination of periodic initial conditions. It shows the intersection point of the orbit in the symmetry plane. Periodic orbits lead to a finite number of intersection points with regular patterns.

Many-body Problem of Gravitation

Chaotic behaviour is not only an interesting theoretical problem, but is of large practical importance, since many phenomena in physics and engineering are described by more than two coupled non-linear differential equations. This includes for example the gravitational processes in three dimensions. In this case the differential equations are non-linear and of the following type for each Cartesian coordinate:

y = -gM y r3 = -gM y x2 + y2 + z2 3 2

Because of the basic type y = -y one expects, that periodic oscillations (orbits) should be possible for certain initial conditions. This is indeed the case for two bodies (in addition there are the cases of a collision for finite size of the bodies and the “scattering” for the case of a body that passes by). For three and more bodies there exist, except for very special initial conditions not long-term periodic orbits, but only more or less chaotic orbits, that can sometimes become quasi-periodic. The apparent regularity of the many-body planet system is a deception. This is due to the relatively short observation time, that is small relative to the time scale, in which the orbits will develop in a chaotic manner.

The situation becomes a bit simpler, if one assumes for the theoretical computation, that all bodies move in one plane, since then the number of coupled differential equations becomes smaller. If one assumes in addition that all bodies have the same size and the same mass m, one can for certain special initial conditions (symmetric configurations) create periodic orbits also for more than two bodies. The following simulation in Fig.9.20 shows such special cases.

With the slider on the left different scenarios can be selected. Thus one can pull individual bodies with the mouse and change the special initial conditions. This leads very soon to a decay of the symmetric configuration. In addition it turns out, that even under such artificial assumptions there exists no long-time stability for more than three bodies, provided the simulation proceeds for a sufficiently long time. The following development can be easily observed, when one zooms into the picture with the slider on the right.

Figure 9.20: Stable and non-stable solutions of the many-body problem of gravitation for the movement in a plane: in the picture an example with five equal masses is shown. With the left slider different initial patterns can be chosen. The right slider allows to zoom in or out of the picture. The initial position of the bodies can be pulled with the mouse.

How could any relatively stable and bound systems develop in the cosmos under these circumstances? One has to consider this as the result of a long-term evolution, with a multitude of collisions and dis-integrations that provide for “friction”, from which those remnants are bound for a limited time to quasi-periodic orbits, that satisfy suitable initial conditions. Most remnants however vanished to the distance, until they interacted again with other systems under the exchange of the energy. On the other hand new candidates entered the original region from other regions leading to a change of initial conditions.

At this point you should again study the essay of Siegfried Grossmann 1 , who studies the question of chaotic systems very throroughly. In the contribution by Guenther Hasinger you can see in much detail, how chaos and collisions can at least lead to order and structure in the cosmos for limited time periods. The simulations by Eugen Butikov simulate a wealth of manybody problems with partially periodic and partially chaotic behaviour.

End of chapter 9