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

Showing changes from revision #7 to #8: 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 let’s look into at some of the software that has been created in theAzimuth code project. There are a number of models there, which are implemented as interactive web pages. The code runs as javascript right in the browser, so their behavior is responsive. This blog covers the stochastic resonance model, by Allan Erskine and Glyn Adgie . The other models also use the same software platform.

### A test drive

The page contains a green sine wave, a randomized, funky curve, and four sliders.

• The sine wave is fed as input to a random process, which generates the other time series.

• Two of the sliders directly control the frequency and amplitude of the sine wave, and affect, in a complex way, the output time series. Verify this through experiment.

• The noise slider controls the randomness of the process. Experiment with it, to get the output to range from perfectly smooth to completely chaotic.

• The “seed” parameter gives different instances of the process.

### Synopsis of the program

This program performs a discrete simulation for a “stochastic differential equation” (SDE), which is an equation that specifies the derivative of a function in terms of time, the current value of the function, and a noise process. When noise is zero, you get an ordinary differential equation.

The program contains the following elements:

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

2. The specific derivative function, which gives the stochastic resonance model

3. Auxiliary functions used by this derivative function

4. Interactive controls to set parameters

5. Plot of the intermediate time series (the sine curve)

6. Plot of the output time series

### Going to the source

I would like everyone now to try to locate the source code, through the following procedure.

• Open the web page for the model.

• Since the code is now running in your browser, you have already downloaded it!

• Now, while your on the page, all you have to do is run your browser’s view-source function. For Firefox on the Mac, use Apple-U, of course, for Firefox on the PC use … (TODO: fill in)

• The window should open up with some html stuff, which after a quick reorganization of your visual parsing neural pathways, should clearly show the text of the web page.

• The code is loaded into the browser by these lines at the head of the html 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>


• Explanation: The first javascript file is for MathJax, which is an open source engine for rendering math formulas. Next, JSXGraph is a cross-platform library for interactive graphics, function plotting, and data visualization. Then we see the main program, StochasticResonanceEuler.js. Finally, normals.js defines a table of random numbers that is used by the calculation.

• These Now four click files on should appear as clickable links. Be warned, however, that the first highlighted two link files, for sadly, StochasticResonance.js look like and they you were should formatted be by there! a robot who is not at the head of his class.

Now that we’re you’re at the source code file forStochasticResonanceEuler.js , I please would humor like me to by ask scanning one through more it, thing of you. Please scan the text of the program, and try attempting to associate sections of the code fragments with elements in the sections of the synopsis that I gave above. You’ll need to put a blur lends over various items of implementation-dependent mumbo jumbo. jumbo; Just just sniff around, and try to form some rough hypotheses.

### Implementation platform

(javascript language, MathJax for formula rendering, and JSXGraph for interactive plotting).

### The stochastic resonance model

Before The proceeding, model works by a simulation that is defined by a stochastic differential equation (SDE), which gives the “instantaneous derivative” of a process as a function of its current value x, time t, and a Brownian-motion-like noise term. Now, to avoid any possible mathematical liability, I must should express point the out following disclaimer. We said that the simulation “infinitesmally is chaotic” driven nature by the “derivative” of the random noise variable term X. raises For challenging ordinary and differential interesting equations questions this for of defining course means the instantaneous rate of change, but which for stochastic equations containing Brownian-motion-like components, the concept involves subtleties that are outside of the present scope scope. of But, this fortunately, article. for But the “Euler-based” discrete numerical approximation, we may get sidestep a these real complications, break here, because the algorithm approximation we are going to use (Euler-based) works by taking a sample of an ostensibly instantaneous rate of change change, whose value may involve random numbers and extrapolates extrapolating that linearly over the interval between sample points. interval; the randomness enters into the picture simply in that the computation of the derivative, at the sample points, may involve random numbers.

Now, In as we said, the derivative stochastic resonance model, the deterministic component of the main derivative variable X is specified as to be a function of time t, the current value x, and a noise term. For stochastic resonance, this Deriv(t,x) is given as the sum of a sinusoidal function of t, called the forcing function, plus a “bistability” function of x, x and that a has random two noise stable variable: points of equilibrium:

Deriv(t, DerivDeterministic(t, x) = SineWave(sine-amplitude, SineWave(t, sine-amplitude, sine-frequency, t) + Bistable(x) Bistable(x), + NoiseSample(noise-amplitude)

where SineWave(t, amp, freq) = … and Bistable(x) =

Here Bistable is the polynomial: Bistable(x) = x * (1 - x^2).

Considered alone, the effect of the sinusoidal forcing function would be to cause the value of X(t) itself to oscillate sinusoidally.

Now what if Bistable was the sole term that defined Deriv(t,x)? Bistable has roots at -1, 0 and 1, so these would be the equilibrium points of X. The equilibria at -1 and 1 are stable, and 0 is an unstable equilibrium which splits the real line into the basin of attraction for -1 and the basin of attraction for +1.

Let’s view each of the basins of attraction as one of the two “states” of a bistable system.

Now consider the effect of the sine function, plus the BiStability polynomial, but still without noise – this is the full deterministic part of Deriv(t, x). If the amplitude of the forcing function is low enough, then, depending on initial conditions, the solution will be drawn towards one of the attractors -1 or +1, and will continue to oscillate around it, as driven by the sine function. So it remains in that state forever.

On the other hand, if the amplitude of the sine wave was large enough, then it would have enough sway to pull the system from one state to the other. The system would then oscillate between the two stable states, at the same frequency as the driving sine wave. It would resonate.

Now, suppose that the sine wave was large enough to periodically pull the system in the neighborhood of the 0, but not large enough for it to actually reach zero and cross over to the other side. It’s still remaining in a single state. What happens, then, if we add some noise to the term for the derivative? Well, if the noise has enough amplitude, and an event comes at the right time, when the system is in the vicinity of 0 due to the sine wave, that event may push the system over to the other side. So we will see random transitions between states, which are most likely to occur at specific phases of the sine wave. With increasing noise amplitude, the transitions will occur with high probability on every cycle, and the system will be resonating with the driving function. The noise here has amplified the effect of the input signal.

### Application to reality: Milankovitch cycles and the frequency of the ice ages

A stochastic resonance model is used in one of the leading explanations for the periods between ice ages.

Here we’ll consider a simplified model, in which the climate of the Earth is modeled as being in one of two possible states: (1) a cold, snowball Earth, or (2) a hot Earth that contains no frozen water. The snowball Earth is in a stable state, because it is white, and hence reflects a lot of solar energy back into space, which keeps it cold. The hot Earth is also stable, because it is dark, which causes absorption of energy, and keeps it hot.

Now we need to add a little bit more realism into the picture, and assume that even in the snowball Earth, the glaciers are concentrated in the northern latitudes, and it is the temperature activity in the northern latitudes that most controls the state of the system. If it gets warm enough up there, the glaciers melt, and the hot Earth state is entered. If it gets cold enough there, then glaciers start to form, and the cold Earth state is entered.

So, we have a bistable system. What about the forcing function?

It turns out that there are astronomical phenomena that cause periodic variation in the amount of solar radiation that is transmitted to the Earth at the northern latitudes. These are known as Milankovitch cycles, and they have periods on the order of tens of thousands of years and upwards. There are three such astronomical cycles: changing of the tilt (obliquity) of the Earth’s axis, precession of the Earth’s axis, and changing of the eccentricity of the Earth’s orbit. (state the separate periods) Now the amplitude of these temperature changes in itself is not sufficient to move the climate from one state to another. (Sound familiar?) But it is possible that random temperature events may cause crossings if they are in sync with the Milankovitch cycles.

Researchers have found interesting correspondences between the actual spacing of the ice ages and the timings of these three Milankovitch cycles.

Well, that pretty much exhausts what I have learned about this subject, so for further information I’ll point you to —-.

### Organization of the program

Now that we have explained the stochastic resonance model, and one of the reasons why it is studied in climate science, let’s return to our regularly scheduled topic, which is the anatomy of this program which implements the model.

The program logic is broken into seven functions.

The top level entry point is the function initCharts, which is short function, and works by dispatching to two other functions, initControls and initSrBoard. (Key: sr = stochastic resonance, board = a graphical gizmo that contains various controls.) As you can see, initControls is building the objects that represent the sliders.

All the logical content of the application is encapsulated in this second function, initSrBoard. Two curve objects are constructed there, called positionCurve and forcingCurve. The forcing curve is constructed from the locally defined forcingFunction, which defines the sine wave. This function reads its values from the amplitude and frequency sliders.

Next, a callback function is attached to the positionCurve object, which makes the call to the stochastic resonance model proper: this is the call to the “mkSrPlot” (short for MakeStochasticResonancePlot) function, which is passed all four of the slider values. This call returns a object that contains a list of time values, and a corresponding list of values for the output variable X.

### The algorithm

The main algorithm is clearly expressed in the function mkSrPlot. Its first step is to construct a function object, called “deriv,” which represents the derivative computation:

Deriv(t,x) = SineCurve + BiStability,

where represents the deterministic component of the derivative.

Then a “stepper” function is constructed, by the call to Euler(deriv, tStep), which returns a function that performs “stochastic Euler extrapolation,” with a step size of tStep, using the given derivative function. The idea of a stepper function is simple: it takes as input the current point (t,x), and a noise sample, and returns the next point (t’,x’). The specific stepper constructed by Euler just increments t by tStep, and increments x by tStep * Deriv(t,x) + the noise sample.

A stepper function of this form is all that is needed for the general toplevel loop, which is implemented by the function sdeLoop, to generate the full time series for the output. This sdeLoop function takes as arguments the stepper function, the noise ampllitude (“dither”), the initial point (t0,x0), a seed value for the randomization, and a number of points to be generated. The loop initializes a currentPoint to (t0,x0), and then repeatedly applies the stepper function to the current point and the next noise sample; the output returned is just this sequence of (t,x) values.

The noise samples are generated by taking a block of values from the (large) array normals(i), and scaling them by the dither value. The “seed” variable controls which section of the array gets used.

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