A non linear equation

Hi guys,

we will consider here a nice non linear equation, which includes (1) DAMPING, (2) HARMONIC OSCILLATIONS, (3) ANHARMONIC OSCILLATIONS, (4) EXTERNAL FORCES. Here is the equation

Before describing the terms in the right hand side one by one let’s spend a couple of words on what is a differential equation.

A differential equation is a relation between (1) the change per unit of time of the state of the system, where the state is simply the pair (position, speed) or (x,v) if you use the symbols shown in the formula, and (2) what the forces in the system are, which is the right hand side. However a differential equation does not have a single solution. We solve an equation like the one above only after we have specified the initial condition. This means that in order to plot the graph of a motion, which shows the dependence (x(t),v(t)) of position and speed upon time, we need to know (1) what is the model, so what is the equation of motion, in our case the one above, and (2) what is the state at the beginning of the motion, therefore (position(0), speed(0)), or (x(0),v(0)). Of course there are properties that we will talk about that are general and hold for every initial state we choose, for example the trajectories of the harmonic oscillator in the phase space are always elliptical etc., but keep in mind that sometimes the graph you are watching can change drastically by changing a bit the initial conditions. That’s the nature of chaos, high sensibility to initial conditions, we will see examples.

Moreover since we will be plotting the motion of the state of the system in time, we will use two kind of graphs. In the first graph typology (time graph) we show the time on the horizontal axis and alternatively the position x(t), or the speed v(t), on the vertical axis. In the second kind of graph (phase graph) we will put the position on the horizontal axis and the speed on the vertical axis. So the motion of the state will be shown as a curve in the plane (x,v) and it is your imagination which produces the motion, because you have no time axis!

To help you imagination I have made this very useful spreadsheet (here is the 97 version) which you can start using straightaway from the next section. The only things I need to tell you are: (1) there are two buttons, “Solve” and “Play” you will be using only “Solve” but check at the instructions to use “Play” at the very end. (2) There are a set of parameters to be used, they are very intuitive and we will not describe them here. Their behaviour will be clear little by little, just follow the instructions below. (3) For the beginning don’t play with the parameter N REFRESH GRAPH, leave it to 10,000!

The Harmonic term

In order to introduce harmonic term ± (ω^2) x, appearing in the right hand side of the equation of motion we will consider a cantilever beam, something which at rest has a vertical position x=0 and which can bend elastically upwards (x>0) and downwards (x<0). If the weight of a 50 tons truck is positioned at the edge of the bridge in the picture, a restoring force equal and opposite is produced “by” the bridge because of the bending. The bending stops only when the restoring stresses exactly balance the weight. The more the weight, the more the bending, so proportionality between the force, which is the acceleration (dv/dt) and bending, which is x

dv/dt = – (ω^2) x

What is going on with the signs? At the beginning we presented the harmonic term as ± (ω^2) x and now only the minus sign is left. The thing is, if you put a + then the restoring force is not restoring anymore, it is diverging and will break the bridge straightaway. What kind of force is it? It is the force you observe when you put a ball on the top of a hill and you gently hit the ball, it will accelerate down the hill’s slope and move faster and faster away. The situation with the bridge is the other way around, here we are talking about the bottom of a valley. If you put a ball on the bottom and you gently displace it, it will experience a restoring force and will oscillate back and forth around the bottom, which is exactly the situation of the bridge. The natural frequency ω is by definition the frequency of this kind of oscillations, provided the oscillations are not too wide. The natural frequency is smaller for elastic materials and higher for hard materials. In terms of sound, the harder the material the more acute is the sound produced by the material when you hit it. The frequency is expressed in Hertz, which is cycle per second. Thus 1 Hz = 1 cycle per s, 440 Hz = 440 cycles per s. Therefore to conclude, when you hit a diapason, you generate vibrations of the solid which have a frequency of 440 Hz (tune A). These vibrations complete therefore 440 oscillations in a second, they are very fast.

To conclude, why did I show the harmonic term with the symbol ±? Because we will consider below the case of anharmonic oscillations where you have (1) a cubic term which effectively restores the motion (has a minus sing) and a competing diverging harmonic term (now plus sign), which dominates for small x. The resulting motion is fascinating, but now let’s stick to the pictures of the motion of the harmonic oscillator. Here is the time graph of position (remember: time t on the horizontal axis and position x(t) on the vertical axis) with initial conditions x(0)=1, v(0)=0, which correspond to displacing the ball (or the bridge) from the equilibrium point by an amount equal to 1 unit of length (x(0)=1) and letting the ball go without giving the ball any speed (v(0)=0)

So, the solution is a cosine function which means that the oscillation goes on forever between x = 1 and x = -1. If you look at the time graph of the speed you will find the same picture shifted by a quarter of a cycle. The starting velocity is 0 rather than 1 as for the position, and the speed oscillates between -1 and 1. Now look at the phase graph, in the following picture a short interval of the motion (the first 5 seconds) is shown (position on the horizontal axis, speed on the vertical axis, the point is moving clockwise)

So both the position and the speed oscillate and the resulting trajectory in the phase space is a circle. That’s because in this example I have set the natural frequency ω to 1, otherwise the resulting trajectory would be describing an ellipse. Here is the phase graph for the same initial conditions but with ω = 1/5 and of course now 25 seconds of motion

If you want to obtain these kind of graphs in motion you can do it with the Excel spreadsheet I mentioned above. All you have to do is to set the parameters to the following values and press “Solve”

Initial point is what we call x(0), whereas initial speed is v(0). Damping is for later, natural frequency is ω, set it to 1 for the circle and to 0.2 to get the ellipse above. Cubic term is for later, ampl ext force and freq ext force for later. Time step is the distance in time of the points on the graph, you can have 1 single jump of 3 seconds or 10 jumps of 0.3 seconds, in the first case do (Time Step = 3, N Steps = 2 (remember that steps are the endpoints of the intervals of motions, so 2 steps make 1 interval!)) whereas in the second case do (Time Step = 0.3, N Steps = 11) you will get the following difference

On the left hand side you have a single point (well two points, if you include the initial point), thus less information about the real trajectory, on the right hand side instead you have 10 points and a smooth trajectory but you need more time to plot them. All the points however lay on the same trajectory, the circle, and the final point is the same for the two trajectories because the same time has elapsed. It is just the interpolation by Excel which produces a line in the first case, there is no line in reality as the figure on the right shows, we are simply closing our eyes and opening them at the end and then tracing a line for convenience, to help the eye.

For those who are interested in the mathematical details, the quality of the solution of the equation is comparable in the two cases. What is the quality? It is roughly the ratio (Time Step)/(Small Time Step), the higher the ratio, the higher the quality. In this case the value of the ratio is 3 / 0.01 = 300 which to me is about 3 times the minimum value one should use.

Now, back to the spreadsheet, as soon as you download it you can open it and press “Solve” because it is already set up for the harmonic oscillator. You can now play with the following parameters: Initial point, Initial speed, Natural Frequency, Time Step, N Steps, Small Time Step. Remember, the ratio (Time Step)/(Small Time Step), at least 100. You should always find nice oscillations, enjoy.

The Damping Term

The damping term - ζ * ω * v is what makes the motion slow down. Suppose that the speed v is greater than zero, then the acceleration (dv/dt) in our general equation of motion (top of the page) has a negative contribution – ζ * ω * v. The contribution is negative because ζ is positive and ω is positive as well. If you are moving with a positive speed and the acceleration is negative then the motion decelerates, you can try and see that the same is true when v is negative. Now, if v is the speed, what are ζ and ω? Let’s start from ω, we know what it is, it’s the natural frequency, why is it there? Forget about it, it is there only to make ζ adimensional, you can consider it just a scaling constant. So concentrate on ζ which in the spreadsheet we call DAMPING. Now, if we add the damping term to the harmonic oscillator equation we get an oscillator with friction, which is a picture much closer to reality

dv/dt = - ζ ω v - (ω^2) x

and you can see that we are getting closer to the final equation. What happens to the graphs of the motion? Take the same values as above and consider three cases for ζ : small, intermediate, large. Here is the time graph we get for the position x(t), with initial conditions x(0)=1, v(0)=0. You can obtain the same graphs in motion if you play with the DAMPING parameter which corresponds to ζ

1) Small ζ (0.03)

2) Intermediate ζ (0.3)

3) Large ζ (1.0)

And here are the corresponding phase graphs

1) Small ζ

2) Intermediate ζ

3) Large ζ

You can clearly see that (1) ζ controls the level of friction in the system, the higher ζ the sooner the system is brought back to the equilibrium (x=0 , v=0); (2) there are two distinct regimes of damping, one (Small ζ) which retains oscillations at a frequency slightly smaller than the natural frequency, and another (Large ζ) where the system is attracted to the equilibrium in one go.

We really talk about an attractor here, what is an attractor? It is a geometrical object, so a collection of points, a circle for example. In our case the attractor is just a single point, the origin (0,0), and the trajectories end up into the attractor because of dissipation.

So, now you can use the Excel spreadsheet and also play with the parameter DAMPING. But now it’s time to add an external force to the system.

External force

Now things really start to become interesting and can get very chaotic. What is an external force? Well, take the case of the tower of a wind turbines, or the mast of your sailing ship if you have one. When the boat is in the harbour and there is absolutely no wind at all then the mast is perfectly vertical. And the response to perturbations is what we discussed above. This means that if you go and manually hit the mast, it will execute damped oscillations. What happens if the ship is moving is, the wind force is steadily applied to the mast and produces a prolonged bending of the mast which steadily counteracts the wind force. This steady bending is the new equilibrium, corresponding to that particular wind speed. If the wind changes the system adjusts with a transient to a new bending, a new equilibrium. Finally If the wind is turbulent we need to add oscillatory terms. A turbulent perturbation is in fact something with eddies, circulating flow of air with oscillatory behaviour.

So to recapitulate, we will look at (1) the case of a constant external force F and (2) of an oscillating external force F sin(Ωt). In order to make things easy we will write the force intensity F as the product (ω^2) X, again, (ω^2) is just a scaling constant, the force is proportional to the dispacement, so remember F = (ω^2) X. It’s just a trick, we control the force via X because in this way we are sure that the number we put in the spreadsheet (X is called Ampl Ext Force) is magically the new equilibrium position that the system will attain under the effect of the force. Let’s see the results then.

Here are the graphs for a constant force X = 10, a damping of 0.1 and initial conditions x(0)=0, v(0)=0

As you can see the trajectory ends up into the new attractor, the point (x=10,v=0). And there is really not much to say a part from the transient oscillations. Like in the case of the mast when the wind changes, because the force is suddenly applied as you can see the system adjusts to the new equilibrium point (x=10) oscillating at (nearly) the natural frequency. Everything else stays the same, try now to change the parameters which we have already introduced together with the new parameter Ampl Ext Force. But remember, in order to have a constant force put Freq Ext Force to zero.

Now, let’s pass to the oscillating force, in this case the expression for the force is F = (ω^2) X sin(Ωt) and we have two parameters to control, the amplitude of the force (Ampl Ext Force, which is X) and the frequency of the force (Freq Ext Force, which is Ω). Now the situation really gets exciting, how do you start playing with the two parameters? Well try, my advise is to forget about the amplitude of the force X, just set it to 1, it only changes the scale of the picture. You should rather play with the frequency. In particular something big happens if the frequency of the driving force coincides with the natural frequency, it’s called the resonance phenomenon. The oscillations get bigger and bigger until the structure breaks down.

You will also notice that if the external force frequency is small the spring gently follows the force, try the following parameters

whereas for high frequency the response is very chaotic, try for example

I will not publish any graph on this topic, try yourself what happens if you change the frequency.

The Cubic Term

The last ingredient of our equation is the Cubic Term, set it to 1 to start with and set the Quad Sign to +1. What this term does is, it changes the picture from a single valley force potential to a double-valley force potential. This produces and interesting behaviour. Start with the following parameters

When you press “Solve” you should observe something like the following graph

What the graph is showing is a ball, which is left free at the point x=2. The bottom of the valley on the right is at the point x=1 so the ball accomplishes a couple of damped oscillations before resting on the bottom valley. If we repeat the same experiment for the point -2 (so change INITIAL POINT to -2 in the spreadsheet) we find that the point is rather attracted to the left bottom valley (x = -1) which is shown in the next picture

So the picture is very simple, two valleys, concurring for which will finally have the ball. Why concurring? Because you have damping, and the solution will be a slowed motion ending pretty soon on a bottom valley. Which bottom valley? It depends on the initial point. You see the picture above, the trajectory describing the ball motion starts at x = -2 and because of the damping it does not succeed in reaching the right valley. But if we give the ball a little help, which means we make the motion start from x = – 2.5 here is the resulting motion

So you see, the ball ends up in the right valley and stays there because its energy is damped out by friction. Nice isn’t it? Well now you can play with the initial condition on the position and found out that the more you put the initial condition far from the origin, the more the number of oscillations before the ball ends up on the bottom valley. Here is what you get setting INITIAL POINT = -6.

Nice and intuitive isn’t it? OK, now let’s go to the juice. The two points x=1 and x=-1 are point attractors for the trajectories of free motion. With free we mean without any external force being applied. Now, when we apply an external force the picture changes completely. There are new attractors. Well for the constant force case the new attractors will be diplaced versions of the point x = 1 and x = -1, just as we get a permanent bending of a cantilever beam when we apply a steady force. You can try yourself, do not exaggerate with the force though, start with little numbers, and of course set FREQ EXT FORCE to zero for a constant force.

Now let’s talk about the nice stuff, the limit cycles. The latter are asymptotical solutions of the motion, thus attractors, which are (1) not made of a single point as before, (2) many in number. Start with the following parameters

as you can see here we start using the N REFRESH GRAPH, it is not 10,000 anymore.. What is that? Since we are going to use Excel in “perpetual motion” we are going to refresh the graph, that is, to delete the old portions of trajectories as the time passes, so that a clean picture of the attractor emerges and not a blackish line out of a cloud of unorganized motion (the initial transient). If you press solve with those parameters you will have a perpetual movie which you will stop pressing ESC and then END. The movie should produce the following cycle (never mind the line in the middle, it’s again Excel’s interpolation), I call it the basic cycle and as you can notice from the list of parameters the amplitude of the force is 0.1, very small, but it’s the frequency which makes possible what we are going to see. I set it to 0.6 so it is neither much greater than 1 (the natural frequency) nor much smaller than 1.

Now, press “Solve”, do you see the cycle forming out of nothing? Cool isn’t it? Now, here is the point, if you increase the force to say 0.11 the limit cycle changes, it increases in size, but qualitatively is identical in shape, have a try. However something topological is going to happen. As we increase the value of the force the limit cycle develops a closed loop, watch out. Here is what you get with AMPL EXT FORCE = 0.13, of course after the transient,

and here is AMPL EXT FORCE = 0.17

For AMPL EXT FORCE = 0.22 the loop is fully developed, although because of the transient we now end up on the left valley, a perfect example of path dependency in chaos theory!

Can you see the loop? I call this diagram the basic cycle with a loop. The transition which is shown in the picture above, which is the transition from basic cycle to basic cycle with a loop is called in Physics a second order transition, something which happens gradually. An example is the behaviour of the residual magnetization versus temperature for a ferromagnet. Everyone should have seen this diagram during the high school, it means that if we leave a magnet in a field and then get rid of the field, the degree of residual magnetisation depends on temperature. How? You can see it yourself int the diagram, it’s very easy, for large temperatures? Say 1000 K, then there is no residual magnetisation, as the story goes entropy wins at high temperatures and order (which means energy minimisation) loses. As we cool the temperature down to about 640 K which is the edge of the curve in the diagram something starts to happen, there is from that temperature downwards always a non-zero residual magnetisation. This is the signal of a phase transition. However it is a gradual transition, you can see it,as a function of the temperature the residual magnetisation doesn’t abruptly pass from 0 to a big value. This is the essence of second order phase transition.

Are we going to explain what a first order phase transition is? Well, of course, because we are about to meet one, let’s push with the gas, which means let’s increase AMPL EXT FORCE, we should be able at some point to see the ball travelling around the two valley shouldn’t we? Possibly it helps to recall here that we have a ball, which is dissipating energy, but we have someone constantly providing energy, an external force, the external force has a magnitude, which is AMPL EXT FORCE and a frequency FREQ EXT FORCE which we have fixed to be 0.6 and we will not touch. So if we increase the magnitude of the force we should be able to sustain a two valleys course for the ball given the friction. Now indeed this is true, try with AMPL EXT FORCE = 0.24, here is what you get

The new limit cycle is what I call the two eyes cycle, it indeed passes along both the valleys. The nice thing about it? The transition between the basic cycle with a loop and the two eyes cycle is a first order transition. What is then a first order transition? Well you might have guessed, it’s the opposite of a second order transition in that in a first order phase transition quantities change abruptly with temperature. Take water and gas, at 99,999 degree C, the system is liquid therefore large density and low entropy, whereas at 100,001 C it is vapor, so low density and high entropy. This is a discontinuity which makes the liquid-gas transition a first order phase transition. In our case, try to plug in these two numbers: (1) 0.235960154208734687619 and (2) 0.235960154208734687620. This is really the limit for Excel, you cannot go beyond, and indeed some computational power might reveal some fractal structure there, but anyway, we will be happy with the two numbers above. What are those two numbers? Well the first produces the basic cycle with a loop whereas the second produces the two eyes cycle. It’s really incredible to watch the movie and see how the first struggles to maintain the two eyes pattern but then dissipates too much and chooses the other attractor, the basic cycle with a loop. And there is no way of going back because too much energy has been consumed.

Wow! So this equation is really hiding things guys! You just have to dig out. For example, we have jumped from AMPL EXT FORCE = 0.17 to AMPL EXT FORCE = 0.22 but of course we have another first order transition there, the one which brings the limit cycle from the right valley, to the left valley, try to find the 21 digit number yourself.

However before leaving you I will present you some final examples. If you move AMPL EXT FORCE up you will see that for certain values you don’t get a limit cycle, but rather something which we would call chaotic. That’s true, if you see something chaotic in life, the first thing you do is to increase the size of the picture, to see if there are stable patterns on a larger scale. You achieve this by increasing the parameter N REFRESH GRAPH. You wil see that indeed, for some values of the force amplitude, patterns emerge and disappear and no limit cycle exist, the attractor is what you see, a changing object, a collection of unstable patterns repeating themselves, you can think probabilistically about this, or you can make an analogy with a resonance in Quantum Mechanics. My message however is that sometimes it’s hard to tell if something is chaos or a transient. Try for example the following parameters, they will also update your N REFRESH GRAPH value, remember, you start the movie pressing “Solve”, you stop the movie with ESC, and then END. If you have the patience to wait until time = 1200 s (on the spreadsheet, 30 seconds in reality) you will see the two eyes cycle emerging out of unstable patterns, give it a go

What sometime happens is that for certain values of AMPL EXT FORCE two competing cycles literally fight with each other, a fight for stability. And we have just seen an example: for AMPL EXT FORCE = 0.45 the two eyes cycle is fighting against other forms of organization and is winning. Let’s see some examples:

1) The fight between the two eyes with a loop and the two eyes with two loops occurring when you set AMPL EXT FORCE = 0.77, (other parameters being the same), at the end you should see the two eyes with a loop shown below emerging

2) The fight between the two eyes with a loop and the two eyes with two loops occurrs for all the spectrum more or less from 0.7 to 1.07. If you set AMPL EXT FORCE = 1.07, you will see that this time the two eyes with two loops is the winner.

3) The edge of chaos at AMPL EXT FORCE = 1.11. It looks like a boiling pot about to explode doesn’t it?

4) The nice multiloop organization at AMPL EXT FORCE = 1.28

Enjoy guys!

PS For what concerns the button “Play”, well, you might want to play around with the other 4 parameters shown on the spreadsheet, they are next to the button “Play”. Keep the value of N STEPS to a maximum of 500 say, because what “Play” actually does is: it (1) calculates the whole trajectory (it overrides N REFRESH and N REFRESH GRAPH) and (2) plots it, (3) waits for a pause (which you specify in seconds) (4) moves the initial position and the initial speed by one step (you can specify the total number of steps of the movie and the size of the steps for both position and speed) (5) go to point (1) again until the end of the steps you have set

…