The Azimuth Project
Blog - The stochastic resonance program (part 1) (Rev #22, changes)

Showing changes from revision #21 to #22: Added | Removed | Changed

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, by Allan Erskine and Glyn Adgie, demonstrates the concept of “stochastic resonance,” which is a widely studied phenomenon that has an application to the theory of ice-age cycles. The Azimuth models are programmed as interactive web pages, which 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 slider controls the frequency of the sine wave, and another controls its amplitude. Try them.

  • The output signal depends, in a complex way, through the “mechanism” of stochastic resonance, on the sine wave. Observe how the amplitude and frequency sliders affect the output signal.

  • The process involves a random component, whose magnitude is controlled by the noise slider. Set it to zero, and see that the output becomes completely smooth. As you increase the noise, verity that the output becomes increasingly chaotic.

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

Inventory of the program

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

Here are the functional components of the program:

  1. Interactive controls to set parameters

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

  3. Plot of the output signal

  4. A function which defines a particular SDE. The stochastic resonance is a property of the solutions to this equation.

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

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.

  • The code is now downloaded and running in your browser.

  • While on 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.

  • Observe 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 load javascript program files, from various locations on the web, into the broswer’s internal javascript interpreter. The first two lines load platform support files: MathJax is an open-source formula rendering engine, and JSXGraph is a cross-platform library for interactive graphics, function plotting, and data visualization. The next line references StochchasticResonanceEuler.js, which is the main code for the model. Finally, normals.js contains a table of random numbers, used in the main program.

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

Now that you’ve reached the source, StochasticResonanceEuler.js, your next challenge is to scan through it and look for associations with the program inventory listed in the preceding section. Try to put a blur lens over any items that look obscure, since the goal here is only to form rough hypotheses about what might be going on.

The concept of stochastic resonance

Now let’s analyze the particular SDE that is being simulated. In this equation, the deterministic part of the derivative is set to a sinusoidal function of time plus a bistable function of the current value:

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

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

Let’a analyze the effects of each of these terms, both separately and in combination.

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

Now let’s consider the bistable polynomial, which has roots at -1, 0 and 1. The root at zero is an unstable equilibrium, and -1 and 1 are points of stable equilibrium. 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 will 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 amplitude is not too large, the system will gravitate towards one of the attractors, and then continue to oscillate around it thereafter – it will never leave the basin. But if it is large enough, the system will be pulled back and forth between the two basins – it will resonate.

Finally, 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 never to actually reach it, or cross over to the other side. If we add in some noise, then a well-timed random event could push the system over to the other side. Therefore, the noise may trigger transitions between states, and this will occur with higher probability at certain phases of the sine wave. Increasing the noise will lead to transitions on more and more cycles of the sine wave will lead to transitions – we may expect to see, in the flip-flopping between states, a “stochastic reflection” of the frequency of the driving sine wave. Even more noise will lead to transitions occurring during a wider range of phases of the sine wave, and a sufficient amount of noise will lead to transitions at all phases – then the noise overtakes the signal, and the output becomes noise.

Moral: under the right conditions, the noise may amplify the effect of the input signal.

Relevance of stochastic resonance to the theory of ice ages

The theory of the timing of the “ice ages” is a fascinating, challenging and unsolved problem in science.

Note: what is colloquial called an ice age – a frozen period – is technically known as a glacial maximum, and the technical meaning of Ice Age is a huge period of time that spans thousands of glacial minima and maxima. There have only been four Ice Ages in the history of the Earth, each of which is characterized by a different pattern in the cycling of the ice ages that it contains.

Here we talk about one of the current hypotheses, which uses a stochastic resonance model. In this account, the climate is modeled as a multistable system, and the forcing is a result of certain cyclical, slowly varying changes in astronomical parameters such as the tilt of the Earth’s axis.

the forcing is produced by astronomical variations known as Milankovitch cycles,

A stochastic resonance simulator program can therefore play a critical role in the testing of such a hypothesis.

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 very simplest model, the climate has two stable – self-reinforcing – 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, and this keeps it hot and melted. The model also includes the fact that the glaciers are concentrated in the northern latitudes, so the northern temperatures are capable of triggering a change in the state of the climate.

There are three astronomical cycles that contribute to the forcing function, all with cycle times measured in in tens of thousands of years:

  • 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 (rotary wobbling) of the Earth’s axis, with a period of 23 kiloyears

Each Together, of these effects produce produces a periodic multi-frequency variation in the amount of solar energy received up in north. the Together, northern they latitudes. sum to a multi-frequency “forcing function.” But the amplitude induced of the temperature changes that are they induce is not large enough to trigger a state change. According to the stochastic resonance hypothesis, it is other, random variations in the state heat change received up north that may be trigger triggered by a random temperature event (of some duration), if it occurs at the right climate phase to change states. One such source of variation is changes in the multi-frequency amount forcing of signal. heat-trapping gases in the atmosphere.

Note that also one such random event could be variations in the amount following of interesting atmospheric interchange greenhouse gases, due to conditions of life on the planet. There is already a consensus that we took have place 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 Azimuth end. blog:

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 Purposes of a stochastic resonance modeling software program

The As we’ve seen, the present program serves to demonstrate the concept of stochastic resonance, and allows you to interactively change explore the its model behavior. parameters and see the results. But what this are kind the broader purposes of such software programs? also has a function within research, both for exploring theoretical questions, and for contributing to the testing of empirical theories.

They For allow you to do a kind case of research, exploring by a means theoretical query in stochastic resonance, suppose we asked how the effectiveness of a forcing function depends on its frequency. This can be explored through a computer experiment, either where you manually vary the frequency parameter, and observe the generated results, or, more systematically, through a meta-program which varies the parameters, and then applies calculated measures to the output signal.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 For starters, we could start replace 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. The model could be used to falsify the theory as well as to validate it.

This is a good case example of the spirit and purpose of scientific programming. programming And at work, and it is an application of Science “science that Really really Matters, 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.

Implementation of the algorithm

Top-level code structure

Our scientific program consists of seven functions. The top-level function is initCharts, which dispatches to initControls and initSrBoard (board means a container for graphical widgets). The job of initControls is to build the sliders.

The program application consists logic of seven functions. The top-level function is initCharts, encompassed within the scope of initSrBoard, which constructs one curve object for the forcing curve, and another for the position curve. Next, the update methods on these objects get set to initControls the appropriate functions. These methods, which are responsible for redrawing the curves whenever its defining input data gets changed, return an object that contains a list of time values and initSrBoard. a (Key: corresponding sr list = stochastic resonance, board = graphical gizmo.) The job of initControls values for the curve. The update method on the forcing curve is set to build a function that computes the sine wave. Note that this function (locally defined) is defined to read the amplitude and frequency values from the sliders. The function “mkSrPlot” gets set as the updater for the position curve object.

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

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 The idea of a stepper function maps is to map the current point (t,x) along and with a noise sample, sample to the next point (t’,x’). The Euler stepper, stepper specifically, maps (t,x) to (t + tStep, x + tStep * Deriv(t,x) + noiseSample).

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

  • The stepper function

  • The noise amplitude, amplitude or (“dither”) “dither”

  • The initial point (t0,x0)

  • The randomization seed offset

  • The Count of the number of points to be generate generated

It The initializes the current point is initialized to (t0,x0), and then repeatedly applies the stepper is repeatedly applies to the current point and the next current noise sample. The output returned is just the 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 The changing randomization offset controls the starting point in the array, we which get leads to different instances of the random process. This is the purpose of the Sample-path variable, which controls the starting point in the array.

Changing the program

How Now that we’ve tried out the program, downloaded its source code, and understood how it works, it’s time to post roll it. up our sleeves and start tweaking it to do new things!

We’re going to proceed by a series of “baby steps.” The first step is to make a local copy of the program on your machine, and get that to run. So, copy the html file and the main java script to a folder on your local machine. I’ll suppose that you’ve stored them into the following folder on your machine: c:\pkg\webmodels.

Now check that the html file is active, by

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