# The Azimuth Project Blog - The stochastic resonance program (part 1) (Rev #46)

This page is a blog article in progress, written by David Tanzer. To see discussions of this article while it was being written, visit the Azimuth Forum. Please remember that blog articles need HTML, not Markdown.

guest post by David Tanzer</i>

At the Azimuth Code Project, we are aiming to produce educational software that is relevant to the Earth sciences and the study of climate. Our programs take the form of interactive web pages, which are written in javascript and use the JSXGraph library for interactive plotting. They allow you to experiment with the parameters of a model and view its outputs. But in order to fully understand the meaning of a program, we need to know about the concepts and theories that informed its creation. So we will also be writing articles to explain the science, the math, and the programming behind these models.

Today I will cover one of the scientific programs developed at Azimuth, written by Allan Erskine and Glyn Adgie. It demonstrates a phenomenon known as stochastic resonance, in which, under certain conditions, a noise signal may amplify the effect of a weak input to a signal detector. The concept was originally applied in a hypothesis about the timing of ice-age cycles, but has since found widespread applications, ranging from the neuronal detection mechanisms of crickets and crayfish to patterns of traffic congestion.

### The mechanism of stochastic resonance

Suppose that we have an input signal which is fed into a signal detection mechanism. The input signal will drive the state of detection mechanism, in a possibly complex, but deterministic fashion. Further, let’s divide the states of the detector into “on” states and “off” states. In the case of a light switch, the input signal may be the force applied to the switch, the state could be the angular position of the switch, and we could say that the sign of the angle determines whether the switch is on or off.

Let’s input a periodic signal, and examine the conditions under which its frequency will be reflected by the flip-flopping of the digital output signal. In particular, let’s consider the case where the input signal is weak, so that the detector remains in a single digital state. How can we amplify such a signal, so that it becomes present in the output? Well, if the input signal is sufficient to drive the detector’s state close to the boundary between digital states, then a bit of random noise may be enough to get it to cross the boundary, into the other digital state. As the input signal moves towards another phase which would have pulled the state further away from the boundary, the system may then cross back into the state where it started. You can imagine that when the settings are tuned just right, the noise will catalyze the detector output to “stochastically oscillate” at the frequency of the input signal.

Stochastic resonance has been found in a wide range of natural and artificial systems, including the signal detection mechanisms of neurons. For instance, there are cells in the tails of crayfish which are tuned to respond to low-frequency signals in the movement of the water, presumably arising from the motions of predators. On their own, these signals are not strong enough to raise the neurons to the firing threshold. But in the presence of the right amount of noise, these signals do cause the neurons to fire.

### Bistable stochastic resonance and Milankovitch theories of ice-age cycles

Stochastic resonance was originally defined in a more restrictive sense, where the two-state system was assumed to be bistable – where each digital state is the basin of attraction for a point of stable equilibrium.

The prototypical application of bistable stochastic resonance was to a form of the Milankovitch theory of the timing of the ice ages.

In the simplest form of Milankovitch theory, the climate is modeled as a bistable system, where one state is a cold, snowball Earth, and the other is a hot, iceless Earth. The snowball Earth is stable because it is white and reflects solar energy, which keeps it cold. The iceless Earth is stable because it is dark and absorbs solar energy, which keeps it hot and melted.

The forcing signals are taken to be the long-duration “Milankovitch” astronomical cycles – tilting of the Earth’s axis, precession of its axis, and variation in the eccentricity of the Earth’s orbit – which vary the amount of solar energy received in the northern latitudes.

The ice-ages are currently running on a roughly 100,000 year schedule, and, intriguingly, the three mentioned Milankovitch cycles have periods of 21, 32 and 100 (FIXME) kiloyears, respectively.

But how would these cycles trigger a change in the state of the Earth’s climate?

In the most recent forms of the Milankovitch hypothesis, a deterministic mechanism is invoked. See, for example Paillard’s paper.

In the original stochastic resonance hypothesis, by Benzi, these Milankovitch cycles produce a signal which in itself is not strong enough to change the state of the climate, but well-timed noise events, at the right phases of a Miloankovitch cycle, could trigger a state change.

The timing of the ice-ages remains a great unsolved problem for modern science.

Stochastic resonance, Azimuth Library

Milankovitch cycle, Azimuth Library

### A mathematical model of stochastic resonance

The demo program uses a stochastic differential equation (SDE) to define a bistable system with a sinusoidal driving function and a noise component. In an SDE, the derivative of the output signal is specified as a function of time, the current signal value, and a noise process.

In this differential equation, the deterministic part of the derivative is set to the a time-varying sine wave plus a bistable function of the current signal value:

DerivDeterministic(t, x) = SineWave(t, amplitude, frequency) + Bistable(x),

where Bistable(x) = x (1 - x2).

Now let’s analyze the effects of these terms.

Alone, the sine wave would cause the output signal to vary sinusoidally.

The bistable polynomial has roots at -1, 0 and 1. The root at 0 is an unstable equilibrium, and -1 and 1 are stable equilibria. The basin of attraction for -1 is all the negative numbers, the basin for 1 is the positive numbers, and the point of unstable equilibrium separates the basins.

We view each basin as one of the states of a bistable system.

Now let’s put the sine wave and the bistable polynomial together. If the wave is relatively weak, the system will gravitate towards one of the attractors, and then continue to oscillate around it thereafter – never leaving the basin. But if the wave is large enough, the system will oscillate between the two basins.

Now let’s consider the noise as well. Suppose the sine wave was large enough to periodically pull the system close to zero – but not enough to cross over. Then a well-timed noise event could push the system over to the other side. There will be a phase-dependent probability of the noise triggering a state change, with more noise leading to transitions on more of the cycles. The flip-flopping between states will contain a “stochastic reflection” of the driving sine wave. Even more noise will cause multiple transitions within each cycle, thereby drowning out the signal and turning the output to noise.

This is the stochastic resonance effect: under the right conditions, the noise may amplify the effect of the input signal.

### Running the demo program

Start by opening the stochastic resonance model web page. On the same plot it shows a sine wave, called the forcing signal, and a chaotic time-series, called the output signal. There are four sliders, which we’ll call Amplitude, Frequency, Noise and Sample-Path.

• The Amplitude and Frequency sliders control the sine wave. Try them out.

• The output signal depends, in a complex way, on the sine wave. Vary Amplitude and Frequency to see how they affect the output signal.

• The amount of randomization involved in the process is controlled by the Noise slider. Verify this.

• Change the Sample-Path slider to get a different instance of the process.

### Going to the source

I would like everyone now to locate the source code, as follows:

• Open the model web page. The code is now running in your browser.

• While there, run your browser’s view-source function. For Firefox on the Mac, click Apple-U. For Firefox on the PC click the right mouse or touchpad button and then select “View Page Source” from the drop-down menu. (TODO: complete for other browsers, Safari, IE, Chrome)

• You should see the html file for the web page itself.

• Note the “script” directives at the head of this file, which refer to javascript programs on the internet. Each line causes the browser to fetch and load a program into its internal javascript interpreter.

   <script src='http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=default'></script>
<script src='http://cdnjs.cloudflare.com/ajax/libs/jsxgraph/0.93/jsxgraphcore.js'></script>
<script src='./StochasticResonanceEuler.js'></script>
<script src='./normals.js'></script>


MathJax is a formula rendering engine, JSXGraph is the graphics library, StochchasticResonanceEuler.js is the model code, and normals.js provides random numbers.

• Click on the link there for StochasticResonanceEuler.js – and you’ve reached the source!

### Anatomy of the program

The program runs a discrete approximation for the stochastic differential equation that was analyzed above. It consists of the following components:

1. Interactive controls to set parameters

2. Plot of the forcing signal

3. Plot of the output signal

4. A function that defines a particular SDE

5. A general SDE simulator, based on the Euler method

The program consists of seven functions. The top-level function is initCharts. It dispatches to initControls, which builds the sliders, and initSrBoard, which builds the curve objects for the forcing function and the output signal (called “position curve” in the program). Each curve object is assigned a function that computes the (x,t) values for the time series, which gets called whenever the input parameters change. The function that is assigned to the forcing curve computes the sine wave, and reads the amplitude and frequency values from the sliders.

The calculation method for the output signal is set to the function mkSrPlot, which performs the simulation. It begins by defining a function for the deterministic part of the derivative:

deriv = Deriv(t,x) = SineCurve(t) + BiStable(x),

Then it constructs a “stepper” function, through the call Euler(deriv, tStep). A stepper function maps the current point (t,x) and a noise sample to the next point (t’,x’). The Euler stepper maps ((t,x), noiseSample) to (t + tStep, x + tStep * Deriv(t,x) + noiseSample).

The simulation loop is then performed by the function sdeLoop, which is given:

• The stepper function

• The noise amplitude (“dither”)

• The initial point (t0,x0)

• A randomization offset

• The number of points to generate

The current point is initialized to (t0,x0), and then the stepper is repeatedly applied to the current point and the current noise sample. The output returned is the sequence of (t,x) values.

The noise samples are normally distributed random numbers stored in an array. They get scaled by the noise amplitude when they are used. The array contains more values than are needed. By changing the starting point in the array, different instances of the process are obtained.

### Making your own version of the program

Now it’s time to start tweaking the program to do new things!

First we’ll make a local copy of the program on your local machine, and get it to run there. Open the html file, and paste its contents into a file, say c:/pkg/stochasticResonance.html. Then paste the source code, and normals.js.

Now point your browser to the file (give file URL example), to make sure that you see the contents, and the model runs. To prove that you’re really executing the local copy, make a minor edit to the html text, and check that it shows up when you reload the page. Then make a minor edit to StochasticResonanceEuler.js, say by changing the label text on the slider from “forcing function” to “forcing signal.”

### Programming exercises

Now let’s get warmed up with some bite-sized programming exercises.

1. Change the color of the sine wave.

2. Change the exponent in the bistable polynomial to values other than 2, to see how this affects the output.

3. Add an integer-valued slider to control this exponent.

4. Modify the program to perform two runs of the process, and show the output signals in different colors.

5. Modify it to perform ten runs, and change the output signal to display the point-wise average of these ten runs.

6. Add an input slider to control the number of runs.

7. Add another plot, which shows the standard deviation of the output signals, at each point in time.

8. Replace the precomputed array of normally distributed random numbers with a run-time computation that uses a random number generator. Use the Sample-Path slider to seed the random number generator.

9. When the sliders are moved, explain the flow of events that causes the recalculation to take place.

I have created the following wiki page for people wish to share any of their code, and notes, related to this article. Just create a section with your name.

### A small research project

What is the impact of the frequency of the forcing signal on its transmission through stochastic resonance?

• Make a hypothesis about the relationship.

• Check your hypothesis by varying the Frequency slider.

• Write a function to measure the strength of the output signal at the forcing frequency. Let sinwave be a discretely sampled sine wave at the forcing frequency, and coswave be a discretely sampled cosine wave. Let sindot = the sum product of sinwave and the output signal, and similarly for cosdot. Then the power measure is sindot2 + cosdot2.

• Modify the program to perform N trials at each frequency over some specified range of frequency, and measure the average power over all the N trials. Plot the power as a function of frequency.

• The above plot required you to fix a wave amplitude and noise level. Choose 5 different noise levels, and plot the five curves in one figure. Choose your noise levels in order to explore the range of qualitative behaviors.

• Produce several versions of this five-curve plot, one for each sine amplitude. Again, choose your amplitudes in order to explore the range of qualitative behaviors.

### Conclusion

Besides serving an educational purpose, such programs also have applications within research: (1) they can be used to perform computational experiments in the pursuit of purely theoretical questions, and (2) they can be used to test empirical theories, by generating their predictions and comparing them to observed and historical data.

For an example of concept exploration, suppose we asked how the effectiveness of a forcing function depends on its frequency. With the current program, we could vary the frequency parameter, and observe the output. More systematically, we could write a driver program that varies the parameters and applies measures to the output signal.

For an example of theory testing, one could imagine a program that implements a model of the Milankovitch astronomical cycles, outputs this signal into the state-changing model of a given theory of climate, and then compares the output signal of the climate model with historical data.

This is the type of scientific programming that we are pursuing at the Azimuth project. It is of interest in itself, being a nice mixture of math, science, and programming, and we will be searching for applications that are relevant to our understanding of the Earth and of life’s problems (what really matters for human survival).

We are starting to prepare for a new round of development at the Azimuth Code Project. Stay tuned, or, better yet, come join us at the Azimuth Forum and help us to work out some specifications for continuation of the work described here.

category: blog