### 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 $Math content$ is $Math content$ 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.

In Fig.9.8 the relative error of the three methods discussed in Fig.9.7 is shown, i.e. for example $Math content$. 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.

#### 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 $Math content$ and for the ordinate the symbol $Math content$.

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 $Math content$. 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.

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

• exponential function $Math content$
• exponential decay $Math content$
• transient processes
• constant velocity $Math content$ with $Math content$ constant.
• constant acceleration $Math content$ with $Math content$ constant.

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

The examples are classified according to the following characteristics:

• divergent (as the exponential function)
• convergent (as the exponential decay)
• periodic
• oscillating and divergent
• oscillating and convergent

In Fig.9.10 a phase space projection $Math content$ versus $Math content$ 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 $Math content$ converges against a finite value, while $Math content$ converges to zero. en

The initial value $Math content$ and the initial abscissa $Math content$ 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 $Math content$ with the differential equation $Math content$. Here the acceleration $Math content$ is increased over eleven steps from $Math content$ to $Math content$. The initial value stays at $Math content$.

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 $Math content$ 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 $Math content$, but this is not the case for $Math content$, 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 $Math content$ at a later time $Math content$. For the numerical solution one then solves the differential equation starting from $Math content$ firstly in the direction of increasing $Math content$ and then in the direction of decreasing $Math content$.

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 $Math content$ and $Math content$ for the start value $Math content$ 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.

$Math content$

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.

$Math content$

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 $Math content$ with $Math content$ and $Math content$ . From the solution for $Math content$ one recovers $Math content$ and $Math content$.

$Math content$

$Math content$

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 “$Math content$”.)

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 $Math content$ shown in the text box the first term $Math content$ is responsible for generating a periodic function and the second term $Math content$ for the exponential decay, as familiar from the differential equations of first order. The factor $Math content$ 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 $Math content$ are added. In the selection box differential equations and initial values for the following functions are preset:

• cosine
• sine
• exponential function
• exponential damping
• hyperbolic sine
• hyperbolic cosine
• delayed oscillation
• accelerated oscillation
• damped oscillation
• growing oscillation

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.

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

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

The three-dimensional phase diagram shows $Math content$ 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 $Math content$, 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 ($Math content$ instead of $Math content$) the differential equation reads: $Math content$.

The term $Math content$ produces a periodic oscillation with period $Math content$, the term $Math content$ an exponential damping and the term $Math content$ a driving force with constant amplitude and the period $Math content$. The very useful step function switches at the given point in time $Math content$ from $Math content$ to $Math content$. 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.

The simulation contains the following pre-defined oscillators:

• free oscillator with adjustable eigen frequency
• dissonant driving force with adjustable frequency
• resonant driving force
• dissonant driving force with damping
• resonant driving force with damping

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

• deflection of a few degrees
• deflection nearly up to the rollover, i.e. angular deflection from the rest point of nearly $Math content$
• shortly after the rollover, i.e. residual velocity at the turning point

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 $Math content$) 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 $Math content$) 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.

#### 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:

$Math content$

$Math content$

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

chaos-theory

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 $Math content$.

##### 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:

$Math content$

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.

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 $Math content$. 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 $Math content$ versus $Math content$ versus $Math content$ 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.

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 $Math content$.

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.

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:

$Math content$

Because of the basic type $Math content$ 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 $Math content$, 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.

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