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

Showing changes from revision #5 to #6: 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 some of the software that has been created in the Azimuth code project , . what it accomplishes, and how it works. There are a number of models there, which are implemented as interactive web pages. The code runs as javascript right in the browser, so they their have behavior a is very responsive. responsive This behavior. blog covers thestochastic resonance model, by Allan Erskine and Jim Stuttard. This will also orient you for looking at the other models, which use the same language and application platform.

Here we’ll focus on the stochastic resonance model, by Allan Erskine and Jim Stuttard. This will also prepare you for looking at the other models, which use the same programming language and application support libraries (javascript language, MathJax for formula rendering, and JSXGraph for interactive plotting).

A test drive

Let’s explore the page. page, There which is contains: a graph showing two functions over time: a green sine wave, and a randomized, funky curve. curve, There are four sliders that control the generated curves. The amplitude and frequency four of the green sine wave are controlled by two of the sliders. These sliders also affect the funky curve. Verify this by experimenting with the sliders.

At the heart of this model is a random process, which generates the second curve. This process is determined by the combination of a deterministic component – known as the forcing function, which here is the sine wave – and a noise component. Experiment with the noise slider, to get the output to range from perfectly smooth to completely chaotic. Try the amplitude and frequency controls, to see how they affect the output signal. Changing the “seed” parameter gives you different instances of the random process.

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

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

  • The noise slider controls how much randomness is injected into the process. Experiment with it, to get the output to range from perfectly smooth to completely chaotic.

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

General Synopsis idea and structure 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. If the noise amplitude is zero, you get an ordinary differential equation, equation. where the derivative depends only on time and the current value of the function.

This The program consists contains of: the following elements:

  1. A general generic SDE simulator (based upon the Euler approximation method)

  2. The formula for the specific SDE used in this stochastic resonance model

  3. Auxiliary functions used by this formula

  4. Interactive controls used to set parameters

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

  6. Plot of the final output time series

Going to the source

I would like everyone now to try to open locate up the source code, and at least skim eyeball through it. Since the code is running in your browser, you have already downloaded it! Here are the steps that I went through to view it:

  • Open the web page for the model.

  • Run your browser’s view-source function. This is browser specific. At home I’m running Firefox on the Mac, and for some reason this function is not appearing in the menu. But a Google search told me to type Apple-U, which worked. (TODO: how to do this with other browsers / environments)

  • Then the window opens up, with some html stuff that clearly contains the text of the web page.

  • But where’s the code? The only possible culprits could be these lines lines: at the head:

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

These lines tell the browser to the following javascript program files:

  1. MathJax.js is an open source engine for rendering math formulas.

  2. JSXGraph is a cross-platform library for interactive graphics, function plotting, and data visualization.

  3. StochasticResonanceEuler.js is the main application code. Click on it!

  4. normals.js, normals.js also clickable, contains an array of random numbers, used in the model

If Now you that had we’re trouble at getting to the file, source you file can get it here StochasticResonanceEuler.js . , I would like to ask one more thing of you. Please scan the text of the program, and try to associate sections of the code with 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, but I know you can do it. Just sniff around, and try to form some rough hypotheses.

Now that we’re there, at the source file StochasticResonanceEuler.js, I would like to ask one more thing of you, before we return to our regularly scheduled programming. Please review the succinct, six-point (precis) structural description of the code, from the previous section. Then skim through this text, and see how many of the elements that I wrote about there you can see being mentioned. This will involve putting a blur lens over all kinds of non-obvious implementation-dependent statements, to sniff out what might be going on there, and form some rough hypotheses.

Implementation platform

The stochastic resonance formula

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

Before proceeding, to avoid any possible mathematical liability, I must express the following disclaimer. We said that the simulation is driven by a function for the “derivative” of the random variable X. Now, for ordinary differential equations derivative this of course means instantaneous time rate of change, but for stochastic equations with Brownian-motion-like components, the definition involves subtleties that are outside of the scope here. But, almost miraculously, for the “Euler-based” discrete numerical approximation that we are using, we are permitted to finagle out of these complications, because the simulator is just a “stepper” that takes an ostensibly instantaneous time rate of change, and extrapolates over the interval between sample points.

The stochastic resonance model

Now, Before as proceeding, we to said, avoid any possible mathematical liability, I must express the derivative following disclaimer. We said that the simulation is driven by the “derivative” of the main random variable X X. is For specified ordinary as differential a equations function this of time course t, means the current instantaneous value rate x, of and change, a but noise for term. For stochastic resonance, equations this containing Deriv(t,x) Brownian-motion-like is components, given as the sum concept involves subtleties that are outside of a sinusoidal function of t, called the forcing present function, scope. a But, “bistability” fortunately, function for the “Euler-based” discrete numerical approximation, we may sidestep these complications, because the algorithm works by taking an ostensibly instantaneous rate of x, change and a whose value may involve random noise numbers variable: and extrapolates linearly over the interval between sample points.

Deriv(t, Now, x) as = we SineWave(sine-amplitude, said, sine-frequency, the t) derivative + of BiStablility(x) the + main NoiseSample(noise-amplitude) variable X is specified as 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, a “bistability” function of x, and a random noise variable:

Considered Deriv(t, alone, x) the = effect SineWave(sine-amplitude, of sine-frequency, the t) sinusoidal + forcing Bistable(x) function + would NoiseSample(noise-amplitude) be to cause the value of X(t) itself to oscillate sinusoidally. What makes the model interesting is the bi-stability term:

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

Let’s Considered see alone, what the effect of the sinusoidal forcing function would happen be if to this cause were the sole value term of that X(t) defined itself Deriv(t,x). to oscillate sinusoidally.

First, Now BiStability what if Bistable was the sole term that defined Deriv(t,x)? Bistable has zeros roots at -1, 0 and 1, so these are would be the equilibrium points of x. X. The derivative equilibria of BiStability is negative at -1 and 1, and positive at zero, so -1 and 1 are stable stable, equilibria, and 0 is unstable. an All unstable equilibrium which splits the points real from line -infinity into to 0 are in the basin of attraction for -1, -1 and all the points from 0 to infinity are in the basin of attraction for 1. +1.

Now, Let’s let’s view consider each the effect of the sine basins function, plus the BiStability polynomial, but still without noise. If the amplitude of the attraction forcing as function one is low enough, then, depending on initial conditions, the solution will converge towards on of the attractors two -1 “states” and 1, and oscillate around it according to the forcing function. If the amplitude of the a forcing bistable function system. is large enough, then the forcing function may have enough activity, at the right moment, to pull the function from one basin of attraction to another.

Next, Now let’s consider the bistability effect polynomial, of plus the forcing sine function, plus the BiStability polynomial, but still without noise – this is the whole full magilla. deterministic Suppose part that of Deriv(t, x). If the amplitude of the forcing function is weak, low and enough, the then, signal is oscillating around an attractor. Then, depending on initial conditions, the noise solution amplitude, a random noise event may pull the signal from one attractor to another. The transitions between basins, then, will contain be both drawn a towards periodic one element, of the attractors -1 or +1, and a will randomized continue element. to oscillate around it, as driven by the sine function. So it remains in that state forever.

This On is interesting / applied in the study other of hand, if the Milkanovich amplitude cycles of the ice sine ages. 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.

Structure of the program

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.

Now the program logic is packaged up in seven functions, and spans a mere 3.5 pages of printed text – it’s impressive that this platform can produce such nice, interactive applications, with this economy of coding.

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

The A top stochastic level resonance entry model point is used in one of the function leading initCharts explanations (hint: for find it and look at it). This is a short little function, which does its work by making two function calls: initControls, and initSrBoard. Please at least scan through initControls, which, as you can see, is building the objects periods that between represent ice the ages. sliders, and returning them as a dictionary of slider objects.

All Here we’ll consider a simplified model, in which the real, climate logical content of the application Earth is encapsulated modeled as being in this one second function, initSrBoard, which is at the end of the file. There, we see that two curve possible objects states: are (1) constructed, a one cold, called snowball positionCurve, Earth, and or the (2) other a called hot forcingCurve. Earth that contains no frozen water. The forcing snowball curve Earth is constructed from the locally defined forcingFunction, which defines the sine wave. Note that in this a function, stable the state, values because of it the is amplitude white, and frequency hence sliders reflects are a used. lot (Note of that solar when energy the back sliders into are space, moved, an event must be fired, which causes keeps the it recalculation cold. to The take hot place. Earth How is this also mechanism stable, implemented because in it the is javascript dark, / which JSXGraph causes application absorption library?) of energy, and keeps it hot.

(How Now we need to make add your a own little version 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 app.) 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.

(Brownian So, motion, we for have noise a process) bistable system. What about the forcing function?

Next, It and turns most out centrally that to there are astronomical phenomena that cause periodic variation in the application, amount a of function solar radiation that is attached transmitted to the updateDataArray Earth member 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 positionCurve tilt object, (obliquity) which of makes the call Earth’s to axis, precession of the stochastic Earth’s resonance axis, model and proper: changing this of is the call eccentricity to of the “mkSrPlot” Earth’s (i.e. orbit. MakeStochasticResonancePlot) (state function, which takes as arguments the values separate for periods) all Now the amplitude of these temperature changes in itself is not sufficient to move the sliders: climate forcing from amplitude, one forcing state frequency, to noise another. amplitude, (Sound and familiar?) seed But value. it This is call possible returns a “plot” object that contains random a temperature list events of may time cause values, crossings and if a they corresponding are list in of sync values with for the main Milankovitch variable cycles. X.

So Researchers now have we found must interesting drill correspondences down between into the logic actual spacing of the mkSrPlot ice function. ages The first step here is to construct a function object, and store in in the variable timings deriv, of which these represents three the Milankovitch derivative cycles. computation:

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 the function Euler(deriv, tStep), which returns a function that performs “stochastic Euler extrapolation,” with a step size of tStep, using the given derivative function function. deriv that was just constructed. The idea of a stepper function is very 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 “dither noise value” ampllitude (amount (“dither”), of noise), the initial point (t0,x0), a seed value for the randomization, and a number of points to be generated. The loop simply 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 the this sequence of (t,x) values values. that are generated (by this iteration).

The noise samples are generated by taking a contiguous block subsequence of values from the (large) array normals[i], normals(i), and scaling them by the dither value. The “seeding” “seed” of the randomization is implemented by using the seed variable as a control value, which controls which subsequence section of the (large) array of normal random variates gets used.

One last point here: How did the Euler function compute a stepper based on what I labeled “stochastic Euler extrapolation?” Well, it’s quite simple: given the deterministic derivative function Deriv(t,x), and step size tStep, it returns the following function:

How to make your own version of the app

((t,x), How noiseSample) to –> post (t it. + tStep, x + tStep * Deriv(t,x) + noiseSample).

(Transitions may resonate with the forcing function, in a stochastic manner)

The time series are generated as sequences of (t,x) pairs. Note that x here is the dependent variable, which is the representation in code of a random variable X.

Problems and challenges

This program could be naturally evolved into a simulator for another SDE, by changing the formula, the intermediate functions, and the associated slider parameters (hint: homework problem to come). Here, the specific components are the SDE derivative formula, the sine wave, which used in that formula, and the sliders for the sine wave parameters.

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

Note also, for a given amplitude of the sine wave, a longer period will have a stronger effect on pulling the state from one pole to another – i.e., the oscillations are more sensitive to low frequencies. One could examine this more systematically by analyzing – or measuring through computational examples – the expected values of the frequency components in the spectrum of the output signal, as a function of forcing amplitude, frequency, and noise amplitude.

  • Effect of frequency

  • How would you study this more systematically, using Fourier analysis in the setting of experimental, stochastic programming

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

  • Modify the code to compute —-

  • 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