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

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>

Today we will look at one of the software models from the Azimuth code project, which aims to produce educational software that is relevant to the study of climate. The program, which was contributed by Allan Erskine and Glyn Adgie, demonstrates the concept of “stochastic resonance.” This is a widely studied phenomenon, which has an application to the theory of ice-age cycles. The Azimuth models are implemented as interactive web pages, which have a responsive behavior because they run right in the browser, as javascript.

### A test drive of the model

Start by opening the stochastic resonance model web page. It displays a sine wave, called the forcing signal, alongside an intricate, time-varying function, called the output signal. There are four sliders, labelled A, B, C and D.

• One of the sliders controls the frequency of the sine wave, and another controls its amplitude. Try them.

• The sine wave is input to a random process, which produces the output signal. Observe how the amplitude and frequency sliders affect, in a complex way, the generated output signal. Can you detect any patterns in the input-output relationship?

• The noise slider controls the randomness of the process. Set it to zero, and verify that the output is completely smooth. See that the output becomes increasingly chaotic as you increase it.

• Changing the Sample-Path parameter gives you different instances of the process.

### Inventory of the program

The program runs a discrete simulation for a “stochastic differential equation” (SDE), which specifies the derivative of a function in terms of time, the current value of the function, and a noise process.

The program contains:

1. A general simulator for SDE’s, based on the Euler method

2. A function which defines the derivative for a particular SDE. The “stochastic resonance” is a description of the behavior of the solutions to this particular equation.

3. Helper functions used by this derivative function

4. Interactive controls to set parameters

5. Plot of the forcing signal (the sine curve)

6. Plot of the output signal

### Going to the source

I would like everyone now to locate the source code, by the following method.

• Open the web page for the model.

• Now, while viewing the page, run your browser’s view-source function. For Firefox on the Mac, it’s Apple-U, for Firefox on the PC it’s … (TODO: fill in)

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

• Note these lines at the top of the file:

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


• These lines refer to javascript program files at various locations on the web. They direct the browser to fetch the code and load it into its internal javascript interpreter. First, MathJax is an open-source formula rendering engine. Second, JSXGraph is a cross-platform library for interactive graphics, function plotting, and data visualization. Third, StochasticResonanceEuler.js is the main code for the model. Fourth, normals.js contains a table of random numbers, used in the calculation.

• Now click on the link for StochasticResonance.js – and you’re there!

Now that you’re at StochasticResonanceEuler.js, your next challenge is to scan through it and try to make loose associations between its terms and clauses and the items of the program inventory given in the preceding section. Try to put a blur lens over any items that look obscure, since the goal here is just to form rough hypotheses about what might be going on.

### The concept of stochastic resonance

In the stochastic resonance formula, the deterministic part of the derivative is set to a sinusoidal function of t, called the forcing function, plus a function of x that has two stable points of equilibrium:

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

where Bistable(x) = x * (1 - x^2).

We consider now the effects of each of these terms, both separately and in combination.

Alone, the sinusoidal forcing function would cause the output signal to vary sinusoidally.

Now let’s analyze the effect of Bistable(x). The roots of Bistable(x) are -1, 0 and 1. This gives stable equilibrium points at -1 and 1, and an unstable equilibrium at 0. The basin of attraction for -1 is all the negative numbers, the basin for 1 is the positive numbers, and the unstable equilibrium point separates the basins.

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

Now let’s consider the combined effect of the sine wave and the bistable polynomial. If the wave amplitude is below a threshold, the system will remain in one basin, forever oscillating around its attractor. But if the amplitude is large enough, it will pull the system back and forth between the two basins – it will resonate.

Now let’s complete the picture by adding in the noise. Suppose the sine wave were large enough to periodically pull the system close to zero, but not enough to cross over to the other side – so it would still remain in a single state. If we add in some noise, then a well-timed random event could push the system over to the other side. The noise, then, will introduce some random transitions between states, with higher probabilities of this occurring at certain phases of the sine wave. As the noise amplitude increases, transitions will occur during more and more cycles of the sine wave will lead to transitions, and the frequency of the flip-flopping will approach that of the sine wave. As noise increases even further, it will dominate the state transitions, which will themselves become noise.

Conclusion: the right amount of noise will amplify the effect of the input signal.

### Relevance of stochastic resonance to the theory of ice-age cycles

A major hypothesis for the timing of the ice ages uses a stochastic resonance model, where the climate is modeled as a multistable system, and the forcing function is produced by astronomical variations known as Milankovich cycles, which periodically vary in the amount of solar energy received at the northern latitudes – where the glaciers are concentrated. The duration of these cycles is on the order of tens of thousands of years, which is at least in the ballpark for the intervals between ice ages.

In the most simplified model, the climate has two stable states: a cold, “snowball” Earth, and a hot, iceless Earth. Each of these states is self-reinforcing. A frozen Earth is white, so it doesn’t absorb much solar energy, and this keeps it cold and frozen. A hot Earth is dark, so it absorbs a lot of solar energy, which keeps it hot and melted.

Now in order to see how the astronomical cycles could be a forcing function that triggers the changes in the state of the Earth, we need to consider something else about the actual Earth: the glaciers are concentrated in the northern latitudes, and an extreme temperature event up north is capable of triggering the climate to change from one state to the other. The snowball Earth is a simplified model for an ice-age climate, and the no-ice Earth for an inter-glacial period. But the stability, or self-reinforcing properties of these models will still apply.

Milankovitch pointed to three astronomical cycles:

• Changing of the eccentricity (ovalness) of the Earth’s orbit, with a period of 100 kiloyears

• Changing of the obliquity (tilt) of the Earth’s axis, with a period of 41 kiloyears

• Precession of the Earth’s axis (rotation of the direction of the axis, like a wobbling top), with a period of 23 kiloyears

Each of these effects produces a periodic variation in the amount of solar energy received up north. Together, they sum to a multi-frequency “forcing function.” But the amplitude of the temperature changes that they induce is not enough to trigger a state change. According to the stochastic resonance hypothesis, the state change may be triggered by a random temperature event (of some duration), if it occurs at the right phase of the multi-frequency forcing signal.

Note that one such random event could be variations in the amount of atmospheric greenhouse gases, due to conditions of life on the planet. There is already a consensus that we have already introduced enough CO2 to stave off the next ice age – but the issue here is we’re maxing it out to the point of disrupting existing life patterns on the hot end.

Now the actual period of the ice ages is on the order of —-, which is in the right ballpark of the periods of the Milankovitch cycles. The question is, how well do the actual specific timings of the ice ages correlate with the multifrequency Milankovitch temperature cycles. This is an open area of scientific investigation.

For further information, see:

### Uses of a stochastic resonance software program

The present program serves to demonstrate the concept of stochastic resonance, and allows you to interactively change the model parameters and see the results. But what are the broader purposes of such software programs?

They allow you to do a kind of research, by means of computer experiment. For instance, one may ask the question: how does the effectiveness of a forcing function depend on its frequency? This can be carried out by means of a computer experiment, either where you manually vary the frequency parameter, and observe the generated results, or more systematically, by means of a meta-program which varies the parameters, and then applies calculated measures to the output signal. This is an example of using computer experiment to address (explore) a theoretical question.

The same modelling software can also be modified to apply empirical tests of theory predictions. We could start in this direction by replacing the sinusoidal forcing function with a function that used computed astronomical data about the Milankovitch cycles. Another step would be to replace the bistable polynomial with a more accurate modelling of the multistable process of the climate. Output could then be “data-mined” for correlations with the actual timings of the ice-ages.

This is a good example of the spirit and purpose of scientific programming. And it an application of Science that Really Matters, which is the motivating principle of the Azimuth Project.

We now return to our study of the particular scientific program at hand, to the anatomy of this stochastic resonance software program.

### Top-level code structure

The program consists of seven functions. The top-level function is initCharts, which to initControls and initSrBoard. (Key: sr = stochastic resonance, board = graphical gizmo.) The job of initControls is to build the sliders.

The overall application logic is contained in this second function, initSrBoard. It constructs two curve objects, forcingCurve and positionCurve. The forcing curve is given by the locally defined forcingFunction, which defines the sine wave. This function reads its values from the amplitude and frequency sliders.

Now a crucial function “mkSrPlot” gets attached to the “update—” method of the positionCurve object. This function is responsible for redrawing the curve, whenever its defining parameters get changed. It returns an object with a list of time values, and a corresponding list of values for X(t).

### The SDE simulator

The simulation is performed by the function mkSrPlot. It first constructs a function for the deterministic part of the derivative:

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

Then a “stepper” function is constructed, by the call to Euler(deriv, tStep). A stepper function maps the current point (t,x) along with a noise sample, to the next point (t’,x’). The Euler stepper, specifically, maps (t,x) to (t + tStep, x + tStep * Deriv(t,x) + noiseSample).

The simulator loop is implemented by sdeLoop, which is passed the following arguments:

• The stepper function

• The noise amplitude, or “dither”

• The initial point (t0,x0)

• The randomization seed

• The number of points to be generated

It initializes the current point to (t0,x0), and then repeatedly applies the stepper to the current point and the next noise sample. The output returned is just this sequence of (t,x) values.

The noise samples are read from an array, called normals, and scaled by the noise amplitude. The array is large, and contains more data points than are needed by the calculation. By changing the starting point in the array, we get different instances of the random process. This is the purpose of the Sample-path variable, which controls the starting point in the array.

How to post it.

### Problems and challenges

• Effect of frequency

• Design a study of the effectiveness of signal transmission, as a function of noise amplitude and signal frequency. How you define the effectiveness measure?

• How would you restructure the code for general, statistical studies of the output time series?

• When the sliders are moved, an event must be fired, which causes the recalculation to take place. How is this mechanism implemented in the javascript / JSXGraph application library?

• Modify to add an exponent slider

• Modify to show graph of expected value (add slider for nTrials) (Not enough random numbers.)

• Add a standard deviation plot

• If you are a climate scientist, let us know of next steps

• Begin to study this book —-, and think of how to write programs for some of the models. Simplify! The hierarchy of models. All models that you post here will be considered as candidates for the Azimuth Code Project page. This may be a way for programmers, ultimately, to give back to the Earth.

category: blog