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

Showing changes from revision #10 to #11: 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 at some software that has been created in the Azimuth 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 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 outputs 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 just 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 theMathJax , which formula is an open source engine for rendering math engine. formulas. Next,JSXGraphJSXGraph* 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.

  • Now click on the highlighted link for StochasticResonance.js – and you should be there!

Now that you’re at the code for StochasticResonanceEuler.js , please your humor challenge me now by is scanning to scan through it, it and attempting attempt to associate make loose associations between code fragments with and elements the points in the program synopsis that given I gave above. You’ll Try need to put a blur lends lens over various any items of that implementation-dependent look mumbo obscure; jumbo; the goal here is just sniff around, and try to form some rough hypotheses. hypotheses about what is going on. heses.

The stochastic resonance model

The In model the works through a simulation defined by a stochastic differential resonance equation model, (SDE), which gives the “instantaneous deterministic derivative” component of the derivative equals a process sinusoidal as a function of its current value x, time t, and called a Brownian-motion-like noise term. Now, to avoid mathematical liability, I should point out that the “infinitesmally forcing chaotic” function, nature plus of the noise term raises challenging and interesting questions for defining the instantaneous rate of change, which are outside of the scope of this article. But we get a real function break here, because the approximation we are going to use (Euler-based) works by taking a sample of an x ostensibly with instantaneous two rate stable points of change, equilibrium: and extrapolating that linearly over the sample interval; the randomness enters into the picture simply in that the computation of the derivative, at the sample points, may involve random numbers.

In the stochastic resonance model, the deterministic component of the derivative is specified to be a sinusoidal function of t, called the forcing function, plus a function of x that has two stable points of equilibrium:

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

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

and Bistable(x) =

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

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

Now let’s analyze the effect of Bistable(x) in isolation. Bistable has roots at -1, 0 and 1, which give the equilibrium points for X(t): -1 and 1 are stable equilibria, and 0 is unstable. The basin of attraction for -1 is all the negative numbers, the basin for +1 is the positive numbers, and the unstable equilibrium is the boundary between the basins.

Let’s treat each basin of attraction as one of the states of the system. It is a bistable system.

Next consider the combined effect of the forcing sine wave and the bistable polynomial. If this the wave amplitude is lower below than a certain threshold, then the system will stay remain in one basin basin, forever, forever oscillating around the equilibrium point. If But the if amplitude it is large enough, it will pull the system back and forth between the two basins, basins i.e., the system will oscillate between the two stable states, at in resonance with the frequency of the forcing function since wave. it will resonate.

Now, suppose the sine wave was large enough to periodically pull the system in the neighborhood of zero, but not enough to cross over to the other side. side So so it still remains in a single state. What If happens if we then add in some noise noise, to then the a term well-timed for random the derivative? Well, if the noise has enough amplitude, and an event comes at the right time, when the system is near zero, that event may push it the system over to the other side. So we will see random transitions between states, which with are higher most probabilities likely to occur at specific certain phases of the sine wave. With As increasing noise amplitude, amplitude increases, more and more cycles of the transitions sine wave will occur lead with to high transitions, probability on every cycle, and the system frequency of the flip-flopping will be approach resonating that with of the driving sine function. wave. The noise here hasamplified the effect of the input signal. As noise increases further, it will become the dominant driver of the state transitions, which will themselves become random noise.

Application Stochastic to resonance and the theory periods of between the frequency of the ice ages

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

In a simplified picture, 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 let’s add a touch more realism to 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 state goes to hot. If it gets cold enough enough, there, then glaciers start to form, and the state goes to cold.

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 (41 kyr), precession of the Earth’s axis (23 kyr), and changing of the eccentricity of the Earth’s orbit (100 kyr). Now the amplitude of these temperature changes in not enough to move the climate from one state to another. (Sound familiar?) But random temperature events may cause crossings if they are in sync with the Milankovitch cycles.

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

That’s as much as I’ve learned about the subject. 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 motivations for studying it, let’s return to our regularly scheduled topic, which is the anatomy of this stochastic resonance program.

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.) Details aside, it is not clear hard to see that initControls is building the objects that represent the sliders.

The main logical content of the application is encapsulated in this second function, initSrBoard. Two curve objects are constructed there, positionCurve and forcingCurve. 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, function with “mkSrPlot” the innocent name of mkSrPlot (for MakeStochasticResonancePlot), gets attached to the “update—” method of the positionCurve object. This function is responsible for redrawing the curve, whenever its defining parameters get changed. The It return returns value an of object this with function contains a list of time values, and the a corresponding list of values output for X(t).

The algorithm

The main algorithm is spelled out in the function mkSrPlot. Its It first step constructs is to construct a function object that gives the deterministic part of the derivative:

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

Then a “stepper” function is constructed, by the call to Euler(deriv, tStep). 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 that is built by Euler maps (t,x) to (t + tStep, x + tStep * Deriv(t,x) + noiseSample).

A stepper function 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. The main loop of the simulator is passed: the stepper function, the noise amplitude (“dither”), the initial point (t0,x0), a seed value for the randomization, and a the 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 taken by from 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 make your own version of the app

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