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

Showing changes from revision #23 to #24: 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 At we will look at one of the software models from the Azimuth code Code project Project , which we aims are aiming to produce educational software that is relevant to the Earth sciences and the study of climate. The Today we will take a look at the stochastic resonance demonstration program, written byAllan Erskine and Glyn Adgie , . demonstrates Stochastic the resonance 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.

Here we will see how to run the program, then explain the idea of stochastic resonance which it demonstrates, and its relevance to a current hypothesis about the ice-age cycles. Then we give a detailed anatomical description of the program, and how it implements the algorithm. We conclude with some challenges for extending the program.

The Azimuth models are interactive web pages. Their behavior is responsive, because they are programmed in javascript, which runs right in the browser. As we shall see, they are programmed using some high-level library support, which makes it easier for the programs to focus on the application logic of the models themselves, rather than upon the presentation mechanisms of the browser. Note that in the future we may decide to scale these programs up from educational demos to larger scale simulations, in which case we will need to develop server-side support for the models.

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 a intricate, chaotic time-varying time-series, function, called theoutput 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 a “mechanism” of stochastic resonance, on the sine wave. Observe Change how the amplitude and frequency sliders to see how they affect the output signal.

  • The amount of randomization involved in the process involves a random component, whose magnitude is controlled by the noise slider. Set it to zero, zero and to see get that a the output becomes completely smooth. smooth As output signal. Verify that as you increase the noise, noise verity slider, that the output becomes increasingly chaotic.

  • Changing Change the Sample-Path parameter gives to get a different instances instance of the random process.

Inventory of the program

The program runs a discrete simulation for a stochastic “stochastic differential equation equation” (SDE), which is specifies a specification for the derivative of the output signal as 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 through the following method. steps.

  • Open the web page for the model. The code is now downloaded and running in your browser!

  • The While code there, is run now downloaded and running in your browser. browser’s view-source function. For Firefox on the Mac, it’s Apple-U, for Firefox on the PC it’s (TODO: fill in)

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

  • The See following lines at the top following of header the lines, html which file load javascript programs, programs from various locations on the web, web into the browser’s internal javascript interpreter:

   <script src=''></script>
    <script src=''></script>
    <script src='./StochasticResonanceEuler.js'></script>
    <script src='./normals.js'></script>

  • Here’s what each line does. MathJax is an open-source formula rendering engine. JSXGraph is a cross-platform library for interactive graphics, function plotting, and data visualization. StochchasticResonanceEuler.js is the main code for the model. And normals.js contains a table of random numbers, used in the main program.

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

Now Your that next you’ve challenge reached is the to source, scan through StochasticResonanceEuler.js StochasticResonance.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, 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 differential SDE equation that is used in the model. This It equation sets the deterministic part of the derivative to a sinusoidal time-varying function sine of wave time plus a bistable function of the current signal value:

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

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

Let’a Let’s analyze the effects of each of these terms, both separately and together.

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

Now 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 stable. 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 stays never in leaves the basin. But if it is large enough, the system will be pulled back and forth between the two basins – it the state will resonate with the driving signal.

Finally, Now, let’s complete the picture by adding in the noise. Suppose the sine wave were was large enough to periodically pull the system close to zero. zero but not enough to cross over to the other basin. If we add in some noise, then a well-timed random event could push the system over to the other side. So the noise may trigger state changes, with higher probability at certain phases of the sine wave. More noise will lead to transitions on more of the cycles – the flip-flopping between states will contain a “stochastic reflection” of the driving sine wave. Further noise will cause transitions across a wider range of phases, and enough noise will drown out the signal, turning the output to 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 open problem in science.

Note: what is colloquially informally called an “ice age” is technically known as a glacial maximum. maximum, The and technical the meaning term of Ice Age is refers to a huge period of time that spans thousands of glacial such minima “ice ages” and maxima. the warm periods between them. There have only been four Ice Ages in the history of the Earth. Each Ice Age is characterized by a different pattern for the glacial cycles that it contains.

One A of the current hypotheses hypothesis for the glacial cycling uses a stochastic resonance model, where the climate is modeled as a multistable system, and the forcing results from certain cyclical, slowly varying changes in astronomical variables such as the tilt of the Earth’s axis. These are known as Milankovitch cycles, and their durations can are be on measured the in scale units of tens 10,000 of to thousands 100,000 of years. years, That which at least puts them in the right same ballpark as the intervals between the “ice ages.”

The Here purpose we of this section is to sketch out the hypothesis – not to make a claim, but rather just to suggest how programs like this can play a role in the scientific enterprise.

In the very simplest model, the climate has two stable states: a cold, “snowball” snowball Earth, and a hot, iceless Earth. Each state is self-reinforcing. A frozen Earth is white, so it doesn’t absorb much solar energy, which keeps it cold and frozen. A hot Earth is dark, so it absorbs a lot of solar energy, and which keeps it hot and melted. The Included in the model also is covers the fact that the glaciers are concentrated in the northern latitudes – that hence thenorthern temperatures are can capable trigger of triggering a change in the state of the climate.

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

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

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

  • Precession (rotary wobbling) of the Earth’s axis, with a period of 23 kiloyears thousand years

Together, These these effects sum to produce a multi-frequency variation in the amount of solar energy received in the northern latitudes. But the induced temperature changes are not large enough to trigger a state change. According to the stochastic resonance hypothesis, it is other, random variations in the heat received up north that may trigger the climate to change states. One such source of variation is changes in the amount of heat-trapping gases in the atmosphere.

Note also the following interesting interchange that took place on the Azimuth blog:

For further information, see:

Purposes of a stochastic modeling program

Our program serves an educational function, which is to show the concept of stochastic resonance, and allow you to interactively explore its behavior. But this type of software also has functions research within applications. research. First, they can be used to empirically explore theoretical questions. Suppose we asked how the effectiveness of a forcing function depends on its frequency. This can be explored, with the current program, by manually varying the frequency parameter, and observing the generated results. On a more systematic basis, we could write meta-program that varies the parameters and applies measures to the output signal.

Such First, software such programs can also be used to generate experimentally explore certain theoretical questions. Suppose we asked how the predictions effectiveness of a theory forcing function depends on its frequency. This can be explored, with the current program, by manually varying the frequency parameter, and compare observing them the to generated actual results. measurement On data. One could imagine, for example, a program more systematic basis, we could write meta-program that implements varies a model of the Milankovitch parameters astronomical cycles, then outputs this signal into the state changing model of a particular theory of climate, and then applies finally measures compares to the output signal signal. of the climate model with observed (or inferred) data.

Such software can also be used to test theories, by generating their predictions and comparing them to actual measurement data. One could imagine, for example, a program that implements a model of the Milankovitch astronomical cycles, then outputs this signal into the state changing model of a particular theory of climate, and then finally compares the output signal of the climate model with observed (or inferred) data.

This is (scientific) programming in the service of our understanding of the Earth.

Implementation of the algorithm

Our scientific program consists of seven functions. The top-level function is initCharts, initCharts. which It dispatches to initControls initControls, which builds the sliders, and initSrBoard initSrBoard, (board which means builds the curve objects for the forcing function and the output signal (called the “position curve” in the program). Each curve object has a container method that is responsible for graphical computing widgets). the (x,t) values for the displayed time series. These calculation methods get called whenever the defining input parameters are changed. The job calculation of method initControls for the forcing curve is set to build a function that computes the sine wave time series. This function reads the amplitude and frequency values from the sliders.

The application calculation logic method is encompassed within the scope of initSrBoard, which constructs one curve object for the forcing output curve, signal and is another for the position curve. Next, the update methods on these objects get set to the appropriate function functions. mkSrPlot, These methods, which are carries responsible out for redrawing the curves actual whenever work its defining input data gets changed, return an object that contains a list of time values and a corresponding list of values for the curve. simulation. The Its update first method action on the forcing curve is set to define 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 deterministic curve part object. of the derivative:

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

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

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

The simulator simulation loop is 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 read from an array normals[i] and scaled by the noise amplitude. The contains many more data points than are needed by the calculation. The randomization offset controls the starting point in the array, which leads to different instances of the random process.

Changing the program

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

We’ll proceed by a series of “baby steps.” First let’s get a local copy of the program to run on your machine. 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