2017-04-29

Chaos game

This post is inspired by the Numberphile video I watched recently: https://www.youtube.com/watch?v=kbKtFN71Lfs.
I greatly recommend watching it and subscribing the channel if you haven't yet done so - it's full of all sorts of mathematical goodness for the mathematical geek in you.

I made an interactive JavaScript implementation of the game, which you can check out HERE (see fig. 1).
The code for the demo is available on GitHub: https://github.com/dagothar/chaos-game. Let me know if you spot any bugs!

Fig. 1. Chaos game demonstration.

So, what is the Chaos Game?
 It's a method for creating fractals. You start with a blank page, a set of vertices arranged in a polygon and a randomly selected starting point. Each turn you select one of the vertices of the polygon and travel towards it covering a certain portion of the distance (for example 1/2 - it is also the same each turn). You draw a dot each time you stop to pick another vertex to travel to. The procedure, simple enough as it sounds, yields some surprisingly complex results.
With a starting set of vertices arranged in a triangle, and covering half the distance to the randomly selected vertex each time, a Sierpiński triangle emerges (see fig. 2).
Fig. 2. Sierpiński triangle is made on a triangular grid with the step parameter set to 0.5.
 More complicated setups yield deservedly rewarding fruits (see fig. 3).
Fig. 3. Pentagonal fractal with step=0.67.


A very curious thing happens when a portion of the game board is barred against the intrusion of the randomly wandering point. In that case, the shape of the restricted area is repeated fractalised all over the board in form of the areas untouched by the point iterations (see fig. 4). The figure shows the blog icon repeated within a square setup of control points. Is it possible to do the same inside a triangle or a pentagon? Does the step size parameter has to be changed?

Fig. 4. DeltaStep logo fractalised in a square after 100,000 iterations.



How to use the demo?
Simply click Start. The game will be played automatically between the vertices of the polygon specified by the red dots. You can Pause the game at any time.
Clear removes the dots, but not the vertices. You can also Download the currently generated image yo your computer.
The Step size slider controls the distance which is covered towards the next randomly selected vertex each iteration (0.0 to 1.0).
You can control the speed at which the new dots are generated. By default, the callback is invoked every 1 ms, but that speed is mostly likely limited by your hardware. You can slow the game down up to 1 dot/sec by using the slider.
You can also customize the colors and the dot size.

Controlling the control points. You can move the control points by hovering over them, clicking and dragging them around. To remove the control point, simply drag it outside the board.
You can add new control points by clicking the Add new point button on the left side of the interface. The new point is placed in the middle of the board, from where you can drag it to wherever you like.
You can make any regular polygon setups by using the Make n-gon button. The number on the right of the button controls the number of the regular polygon vertices.

Defining the restricted area. You can paint on the board with your mouse. Simply click and drag around outside the control points. The painted area will be blocked for the game iterations, thus allowing you to create the fractals of custom patterns. (Note: do not block too much of the board - the script may get stuck trying to find a new possible jump target).

There are several unanswered questions the game leaves us with. What happens when the control points are placed on the vertices of a concave polygon? What if the polygon is distorted? What does the step size actually do? How do you achieve the fractalisation effect with more or less points than four?

The chaos game invites.


Literature
[1] Numberphile, Chaos game video, available: https://www.youtube.com/watch?v=kbKtFN71Lfs
[2] Wikipedia,  Chaos game, available online: https://en.wikipedia.org/wiki/Chaos_game

2017-04-23

Dragon curve

Fig. 1. A forest...? a whirlpool...? a blaze? of Dragons. What do you think a collective noun for these should be?
You fold a paper strip absent-mindedly and a wondrous creature emerges...

If you read The Jurassic Park by M. Crichton (and I’m sure you did), you sure have noticed how each of the chapters started with a page featuring an interesting drawing. It began at first with a simple set of lines, but as you progressed in the book, it got more and more convoluted. I couldn’t figure out the rule behind it at the time, and I found the explanation only years later.

Fig. 2. Dragon curve animation by 碳酸鈣 (available at  this
 Wiki page under CC license).
The lines form a fractal, which is also called Heighway dragon, the Jurassic Park dragon, or simply a dragon curve (see fig. 2). It’s quite a nice thing to doodle! There are a few ways you can do this and the choice is yours:


1. The paper strip folding method
Take a strip of paper. Start folding it: right over left, and keep repeating: right over left, right over left… You won’t be able to do that more than 4-5 times. Don't worry! You can always make more components and connect them end to end  to get more complicated shapes (more iterations of the fractal).
Unfold  each of the components and crease all the folds so they all form right angles (see fig. 3). Place it on a flat surface, and voila!

Fig. 3. How to make the Dragon Curve by folding the paper strips.


2. The pen-and-paper method
You can very well draw the dragon on the paper, starting with a straight line segment and making proper right or left 90 degree turns. There is a neat method that can be used for the determination of the turns in the sequence. Let's start with an example (R marks a right turn and L - a left turn):

1st iteration:  R,
2nd iteration:  R_R_L,
3rd:            RRL_R_RLL,
4th:            RRLRRLL_R_RRLLRLL,


And so on… As you can see, in each of the iterations you take the previous one, add an R (right 90 degree turn) at the end, and then add the previous iteration again, but this time  both reversing i, and swapping R’s and L’s (see fig. 4).
Fig. 4. Drawing the Dragon Curve.


3. The computer friendly method
You can calculate the n-th turn direction by evaluating the following expression:

if ((n & -n) << 1) & n == 0 then:
  turn ← R;
else:
  turn L;

& is a bitwise and operator, << is left shift and you have to remember the -n should be encoded as two-complement. Also, the n-th full iteration is a sequence of 2n-1 turns, so you should probably calculate all of these! :)

For example, the direction of 42nd turn is:

((00101010 & 11010110) << 1) & 00101010 =
= 00000100 & 00101010 = 0,
so the 42nd turn is R.

The figure 1 comes from an old program I wrote in C to explore the possibilities offered by the Dragon curves. You can find this program here.

If you ever want to summon a dragon, just use your imagination



Note. This was originally posted on my Tumblr blog: http://dagothar.tumblr.com/post/26697353061/the-jurassic-park-fractal-if-you-read-the. I think this is a much better place for it though, and I'm pretty sure you haven't seen it yet anyway!

2017-04-21

Quick addendum: Lotka-Volterra equations

In my previous post about the simulation of a simple virtual prey-predator  ecosystem on the world of Wa-Tor, I have mentioned Lotka-Volterra equations as another simple model for such relations. I think it would be interesting to observe how this model performs - is there something curious these equations can show us?

Just to remind - Lotka-Volterra equations describe a relationship between the sizes of two populations: prey (denoted as x) and predators (denoted as y). The equations are:

ẋ  = αx - βxy 
ẏ  = δxy - γy

The rate of change of the prey population () depends on the natural rate of reproduction (αx) and the rate of killing by predators (-βxy). The latter is modelled based on the chance that the animals of two species encounter each other - the more of the either, the higher the chance! The predators benefit from those encounters, and thus gain numbers at the rate described by the term δxy. The predators, being at the top of the food chain can only die to accidents, sickness and the old age. This is modelled by the term proportional to the current predator population (-γy).

In this interactive chart below you can check the various outcomes for such an ecosystem by varying the parameters of the equation (α, β, γ, δ, and the initial populations x0 and y0):


L-V parameters Initial conditions sim. parameters

The Lotka-Volterra model, besides having quite a bad-ass name, is quite powerful. Still, it suffers from the problems caused by its simplicity. One of them is known as the atto-fox problem. Notice how the population numbers do not actually describe the size of population in the units of individuals, but rather in terms of relative sizes. Yet, it often happens that due to a choice of parameter, one of the population sizes (say that of the predators) dwindles to obscenely low values - there would be less than one animal left! As an example, only 10−18 of a fox may have been left alive at some point in time. Obviously, any real population cannot bounce back from such a disaster. Thus the model reaches its boundaries.

2017-04-19

Wa-Tor

(Check out the simulator HERE)

Fig. 1. Wa-Tor simulator.
Introduction
I could not give an introduction to the world of Wa-Tor more perfectly than by quoting its inventor, Alexander Keewatin Dewdney [1]:
Somewhere, in a direction that can only be called recreational at a distance limited only by one's programming prowess, the planet Wa-Tor swims among the stars. It is shaped like a torus, or doughnut, and is entirely covered with water. The two dominant denizens of Wa-Tor are sharks and fish, so called because these are the terrestial creatures they most closely resemble. The sharks of Wa-Tor eat the fish and the fish of Wa-Tor seem always to be in plentiful supply.
Wa-Tor is a simple model of population dynamics in the context of two interacting species: one of them is prey and the other are predators (in Dewdney's implementation they are respectively: fish and sharks).

We know how that goes: in the time of abundance, the predators multiply without second thought - so much that the food now becomes depleted and they die of starvation. The prey can now graze undisturbed and returns to its high numbers quickly - and thus the cycle repeats.

A simple model for that process is given by a set of differential equations: Lotka-Volterra equations. The equations are as follows:

dx/dt = αx - βxy
dy/dt = δxy - γy

The x represents the number of prey. In time, it increases due to breeding and the number of current prey (thus the term +αx). The prey is eaten by the hunters, and the rate is dependent on how often the two species meet (the term -βxy). Similarly, the numbers of predators y grow due to the food supply (+δxy). Hunters die though - either due to accidents, or due to starvation, which is refleted by the last term: -γy. A sample solution to Lotka-Volterra equations is presented in fig. 2 (the Matlab code is also available in the simulator repository).

Fig. 2. Lotka-Volterra solution for α=2/3, β=4/3,
γ=1 and δ=1.
 
While the ecosystem described by Lotka-Volterra model never  crashes, the Wa-Tor presents a different story. It is more than likely that a selection of its parameters will result in one -or both- populations dying out. It is quite tricky to find the values that make the ecosystem of the planet stable (I already put the correct values in the simulation; if you feel adventurous you may change them).

Conveniently, the original paper of Dewdney [1] is available online, and it's sure worth taking a look at!


Rules of Wa-Tor
I shall give a brief presentation of the rules of Wa-Tor here. They are described in much finer detail in Dewdney's paper [1], but stay away from Wikipedia article (it has some bugs)!

The planet of Wa-Tor is represented as a two dimensional array of cells (squares). The array has no boundaries: if you pass through the left edge you wind up on the far right side (and vice versa). The same applies to the top and bottom edges. While boundless, the planet has only a limited surface, measured by its width and height.

The time of the planet Wa-Tor is discrete: it passes in turns (called chronons by Dewdney - it's such a cool word!). Each chronon the state of the world is updated:

  1. Each fish (prey) moves. The fish move randomly, and only to an empty adjacent square. If no free space is available, the fish stays put. If the fish had reached its reproduction age while moving, it leaves offspring (a new fish) in the space it has just left. The age is then set to 0.
  2. Each shark (predator) moves. If there are any adjacent fish, the shark picks one at random and devours it. If there are no fish in the neighbourhood, the shark moves randomly, just like the fish. If the shark has reached its reproduction age, it leaves a new shark in the square it left behind. The age is then set to 0. Each turn the shark gains an unit of starvation. If a certain starvation level is reached, the shark dies. If a shark eats a fish, its starvation level is set to 0.
  3. The screen is updated!
These steps are repeated over and over.
When speaking of adjacency, the classic Wa-Tor algorithm uses the von Neumann idea of the neighbourhood: only the squares that are directly left, right, up and down from a given square count as neighbours.

The parameters are then as follows:
  • width - the size of the world in x direction,
  • height - the size of the world in y direction,
  • initial prey - the number of fish at the start of the simulation,
  • initial predators - the number of sharks at the start of the simulation,
  • prey reproduction age - the number of turns after which a fish produces an offspring,
  • predator reproduction age - the number of turns after which a shark produces an offspring,
  • predator starvation age - the number of turns without finding prey after which a shark dies of hunger,
  •  optionally: age variance - a parameter that controls how random the initial ages and starvation levels are for newly born animals.


Wa-Tor Simulator
I made a simulator of the Wa-Tor world in JavaScript (see fig. 1). You can access it HERE. The code is available on GitHub: https://github.com/dagothar/wator.

I hope I made the simulator robust enough and parametrizable enough! You can probably change any relevant parameter (I skipped on the color customization this time though. I also wanted to put in sprites to represent the species but decided to abandon the idea. They wouldn't be readable enough on larger boards.)

Here's a list of features:
  • Wa-Tor simulation according to the original rules (plus a few modifications),
  • parametrizable: set the world size, age & starvation parameters for both prey and predators, initial species distribution,
  • watch the real-time data on the number of the prey and the predators,
  • a dynamically updated chart which shows the population stats,
  • optional: imagine the blue squares are people and the red ones are raptors!
Note: You have to 'Reset' every time you change the parameters.

It is quite captivating to observe the ecosystem evolve in time. The fish gather in groups which are then dispersed by sharks.

Here is how it looks like!
Fig. 3. An era in the world of Wa-Tor.
Fig. 4. This is how the population of the species change. The scale for the fish is on the left, the scale for the sharks - on the right.
An insightful observer may perhaps find the dynamic world of Wa-Tor a valid parallel for realities of human existence.

Notes on Implementation
The inventor of Wa-Tor supplies an abundance of pointers in his paper on how to implement the model. I tried to keep very close to the original idea, but included some changes of my own that I think improve on the performance and aesthetics of the simulation.

Dewdney suggests using one of two data structures for holding the information on fish and sharks: either a set of tables, or a set of linked lists. Both have their advantages. I use the lists to iterate through the animal instances in the update step, while still keep the tables to make the search for adjacent empty squares/prey faster.

The original algorithm calls for iterating through the table of animals each chronon and updating their state/position in turn: from the top left to the bottom right corner. This is not ideal, since it introduces artifacts in how the animals can move and group together - the ones that move first are of an advantage. I shuffle both the prey and predator lists before making the update step, so that the movement seems much more natural.

It also looks more natural when the update step for the predators is done before the prey moves. In this way, you can see the predators actually moving to the squares occupied by prey. If it were the other way around, the predators would seem to move to the squares occupied by updated prey and which were not yet displayed on the screen.

The author does recognize that while seeding the Wa-Tor world with the initial distribution of animals it is best to slightly randomize their age, so that once they reach maturity the population does not suddenly double in size. I do that by introducing a random age offset when placing new animals (also when they are born!). The offset distribution is centered at 0, so that no bias is introduced. I also vary the starvation age of the predators in the similar fashion.

Additionally, I made it possible to use the other classic rule of the neighbourhood for the model: Moore neighbourhood (4 squares on the: left, right, top, bottom + diagonals).



Can you play God and set the rules of the planet Wa-Tor such that it thrives for eternity?
Your chronon!


Literature
[1] A. K. Dewdney, Sharks and fish wage an ecological war on the toroidal planet Wa-Tor, Scientific American, December 1984, pp. 14-22, available online: http://home.cc.gatech.edu/biocs1/uploads/2/wator_dewdney.pdf

2017-04-16

Computus

Easter is the most important Christian holiday which celebrates the the resurrection of Christ from the dead (strangely enough, the name for the festival comes from the name of the long forgotten Germanic goddess of the radiant dawn: Ēostre). It is a movable feast - and the task of determining the precise date caused much problems for the monks and scholars over the ages.

The task was called Computus - from latin: the calculation, and it was indeed the calculation of its age. The story of how the Easter date was computed is quite curious and I recommend reading the very interesting article on the Wiki: link.

One of the algorithms for the Easter date computation was proposed by the famous Carl Friedrich Gauss in 1800. The algorithm is:

function easter(year):
  a ← year mod 19
  b ← year mod 4
  c ← year mod 7
  k ← floor (year/100)
  p ← floor (13 + 8k)/25
  q ← floor (k/4)
  M ← (15 − p + k − q) mod 30
  N ← (4 + k − q) mod 7
  d ← (19a + M) mod 30
  e ← (2b + 4c + 6d + N) mod 7

  if 22+d+e <= 31:
    month ← 3
    day ← 22+d+e
  else:
    month ← 4
    day ← d + e - 9
    if day == 26:
      day ← 19
    end
    if day == 25 and (11M + 11) mod 30 < 19:
      day ← 18
    end
  end

  return (year, month, day)

It should be easy to code that down in any language you know. Sample JS implementation is available below. Isn't it nice that the next year's Easter falls on the Fools Day?

Date of the easter for the year is: (date here)


HAPPY EASTER!

2017-04-14

Lissajous

Lissajous curves are a family of curves which can describe a complex harmonic motion. They are described with the following parametric equations:

x(t) = Ax sin (ωx t)
y(t) = Asin (ωt + φ)

The parameters Ai define the amplitudes of the oscillations, the ωi are the frequencies and the φ is the phase offset between the two. The curve is typically ploted through the whole rotation: t = 0, ..., 2π.

The Lissajous curves could be produced by harmonographs (an apparatus that has a pendulum with some writing device attached to it), or on oscilloscopes, where the horizontal and vertical electromagnets would drag the electron beam along the screen. The general shape of the curve is very much dependent on the ratio between the frequencies ωx and ωy. The ratio defines the time it takes for the two orthogonal oscillations to coincide and thus the curve itself becomes more crossed and twisted. The change in phases φ in turn makes the curve rotate.

Lissajous curves are typically defined in two dimensions, but nothing bars an extension into 3D. In three dimensions the equations become:

x(t) = Asin (ωt + φx)
y(t) = Asin (ωt + φy)
z(t) = Asin (ωt + φz)

The parameters are quite alike as the ones described above. The only difference is the introduction of separate phase offsets for different axes.

I was quite interested to see the Lissajous curves in 3D. The changing of phase parameters while considering 2D equations often results in the plots twisting as if they were only a projection of a higher-dimensional entity which rotates. Adding the third dimension we can confirm it is indeed the case.

I made a Lissajous 3D curve generator, which is available HERE
The source code can be found on GitHub: https://github.com/dagothar/lissajous

The application is shown in fig. 1 below. You can change the parameters of the curve on the right: the amplitudes and the frequencies for the respective axes (X, Y and Z) can be selected with the combo boxes, while the phase shifts are set using the sliders. You can rotate the 3D view on the left by dragging the mouse around it.

Fig. 1. Lissajous curve generator.


And here's a small gallery of examples.

Fig. 2. The most basic curve is just a circle (ωx=1, ωy=1, ωz=1).
Changing the phase turns it into an ellipse.
Fig. 3. The knot seems to rotate in the plane when changing the phase parameter  (ωx=2, ωy=3, ωz=0, φx=0, φy=0:2πφz=3.14).
Fig. 4. The same knot as above with the third dimension added. The 3D curve projects into the same image.
The rotation is now clearly visible.
Fig. 5. Changing the ratio between the frequencies ωx and ωy.

Fig. 6. This one looks like a twister (ωx=9, ωy=10, ωz=1, φx=0, φy=1.57, φz=3.14).

Fig. 7. Orbits? (ωx=24, ωy=25, ωz=22, φx=0, φy=1.57, φz=3.14).

Some of these shapes are quite surprising. What interesting and beautiful shapes can you find?

2017-04-12

Numerical integration

Differential equations are amazing and useful stuff. Not a day passes by when you don’t solve at least one. And if you do not - don’t you feel like you have wasted a day of your life?
Sadly, as much as you’d wish to, you don’t always have your trusty Matlab (or Mathematica etc.) with you. Never fear, a sheet of paper (and later - a spreadsheet) will be more than enough!

The Euler method

You already know how to solve differential equations in memory:
A car is travelling with the velocity of 10 m/s and is now at the 100 m mark on the race track. What will its position be after 1.5 seconds?
The answer is of course: 115 m.
What took the derivative of the position of the car (the velocity) and added it to its current position after multiplying it by the time step (1.5 s). This is the essence of the Euler’s method of numerical integration:
x_{n+1} \; = \; x_{n} \: + \:  \Delta t \: f \: (t_n,\: x_n) \\
t_{n+1} \; = \; t_{n} \: + \: \Delta t
Basically, you add to the current state (the “position”) the derivative of the state (the “velocity”) multiplied by the length of the time step. If the steps are small, the result should be good enough.
Let’s solve a simple differential equation numerically (you might recognize this as an inertial element):
\dot{x} \; = \; u \: - \: \frac{1}{T} \: x
y \; = \; \frac{k}{T} \: x
We do it in the following steps:
1. Prepare the table. It should have a column for time (t), the input (u), the state (x), the state derivative column (dx) and the output column. Fill in the time values at desired time steps.
t u x dx y
0.0
0.1
0.2
2. Next, fill in the input values at different values of time. Here, we can use an unit step for example (so all the input values will be 1). We can also pick the initial value for the state. Let’s say it’s 0.
t u x dx y
0.0 1 0
0.1 1
0.2 1
3. Let’s start. First, we calculate the derivative dx at t=0.0. We use the equation \dot{x} \; = \; u \: - \: \frac{1}{T} \: x inserting the values from the current row: 1\:-\:0\;=\;1. We can calculate the output at the same time: y\;=\;x\;=\;0.
t u x dx y
0.0 1 0 1 0
0.1 1
0.2 1
4. For the second row, we have to update the state. We add the dx from the row above, multiplied by the step \Delta t: 0 \:+\:0.1 \: \cdot \: 1\;=\;0.1.
t u x dx y
0.0 1 0 1 0
0.1 1 0.1
0.2 1
5. We repeat the step 3 for the second row now, calculating new dx and the new y.
t u x dx y
0.0 1 0 1 0
0.1 1 0.1 0.9 0.1
0.2 1
6. Keep repeating steps 3 and 4 for the subsequent rows:
t u x dx y
0.0 1 0 1 0
0.1 1 0.1 0.9 0.1
0.2 1 0.19 0.81 0.19
0.3 1 0.27 0.73 0.27
0.4 1 0.34 0.66 0.34
0.5 1 0.41 0.59 0.41
0.6 1 0.47 0.53 0.47
0.7
(If you make a plot of this it should look like the one in fig. 4.)
Solving various differential equations make for an excellent pastime for the boring and useless meetings and seminars. You will rarely get scorned for doing maths on such occasions; certainly not as if you were doing sudokus, crosswords or mind-absently doodling.
Solving differential equations numerically by hand will command great respect among your peers. Whenever I explained what I was doing to a curious inquiring soul, they never bothered me with questions again.
They evidently understand how important math is!

The Runge-Kutta method

The Euler method is fast and easy to use - and certainly sufficient for simple demonstrations. It is not very precise though.
Another popular method that helps with that problem is the classic Runge-Kutta fourth-order method. In RK4, the next state value is calculated as:
x_{n+1} \; = \; x_{n} \: + \: \frac{\Delta t}{6} \: (k_1 \: + \: 2k_2 \: + \: 2k_3 \: + \: k_4) \\
t_{n+1} \; = \; t_{n} \: + \: \Delta t
k_1 \; = \; f \: (t_n,\: x_n)
k_2 \; = \; f \: (t \: + \: \frac{\Delta t}{2}, \: x_n \: + \: \frac{\Delta t}{2} \: k_1)
k_3 \; = \; f \: (t \: + \: \frac{\Delta t}{2}, \: x_n \: + \: \frac{\Delta t}{2} \: k_2)
k_4 \; = \; f \: (t \: + \: \Delta t, \: x_n \: + \: \Delta t \: k_3)
We simply calculate the derivatives four times. First - based on the current time and state values. The next approximation is based on a what-if scenario: what if we used thus calculated state value to compute derivatives at half the time step from now (t\:+\:\frac{\Delta t}{2})?
The third coefficient is calculated the same way, but taking the newly made prediction into account. At last, the what-if derivative is calculated using k3 prediction at the end of the time step (t\:+\:\Delta t).
A simple implementation of RK4 in JS might look like this:
function solve(fun, t, u, x, dt) {
      // calculate derivatives k1
      var k1 = fun(t, u, x);

      // update state for calculating x2
      var x1 = x.slice();
      for (var i = 0; i < x.length; ++i) {
        x1[i] += k1[i] * dt/2;
      }

      // calculate derivatives k2
      var k2 = fun(t+dt/2, u, x1);

      // update state for calculating x3
      var x2 = x.slice();
      for (var i = 0; i < x.length; ++i) {
        x2[i] += k2[i] * dt/2;
      }

      // calculate derivatives k2
      var k3 = fun(t+dt/2, u, x2);

      // update state for calculating x4
      var x3 = x.slice();
      for (var i = 0; i < x.length; ++i) {
        x3[i] += k3[i] * dt;
      }

      // calculate derivatives k4
      var k4 = fun(t+dt, u, x3);

      // update state with the final value
      var nx = x.slice();
      for (var i = 0; i < x.length; ++i) {
        nx[i] += dt/6 * (k1[i] + 2*k2[i] + 2*k3[i] + k4[i]);
      }

      return nx;
};
You could use it like this:
// a function that calculates the derivatives
function odefun(t, u, x) {
    var dx = [ 0.0, 0.0 ];
    dx[0] = x[1];
    dx[1] = u - 0.1*x[0] - 0.1*x[1]; 

    return dx; 
};

var x = [ 0.0, 0.0 ]; // initial conditions
for (var t = 0.0; t < 10.0; t += 0.1) {
    u = 1.0; // the input
    x = solve(odefun, t, u, x, 0.1);
    console.log(t, x);
};

RK4 spreadsheet

While Euler’s method is perfectly portable and easy on the hand while drawing it out on the paper, RK4 is somewhat more complicated. I made a little spreadsheet that can be used to solve simple LTI systems in LibreOffice Calc or Excel.
The spreadsheet (and other materials mentioned here) is available HERE
LTI systems (linear time-invariant systems) are a class of systems taking the following form:
\mathbf{\dot{x}} \; = \; A \: \mathbf{x}(t) \: + \: B \: \mathbf{u}(t)}  \\ \mathbf{y} \; = \; C \: \mathbf{x}(t) \: + \: D \: \mathbf{u}(t)
A, B, C, D are matrices here, which describe the linear relationships between the input signals u, states x and output signals y (x, u and y are vectors). The above equations are presented in condensed notations, but are really quite simple. To be more verbose:
\dot{x}_1 \; = \; A_{11}x_1 \: + \: A_{12}x_2 \: + \: ... \: + \: B_{11}u_1 \: + B_{12}u_2 \: + \: ...
\dot{x}_2 \; = \; A_{21}x_1 \: + \: A_{22}x_2 \: + \: ... \: + \: B_{21}u_1 \: + B_{22}u_2 \: + \: ...
...
y_1 \; = \; C_{11}x_1 \: + \: C_{12}x_2 \: + \: ... \: + \: D_{11}u_1 \: + D_{12}u_2 \: + \: ...
y_2 \; = \; C_{21}x_1 \: + \: C_{22}x_2 \: + \: ... \: + \: D_{21}u_1 \: + D_{22}u_2 \: + \: ...
The A matrix (state matrix) describes how the states influence the state derivatives. The B matrix (input matrix) relates the influence of input signals on the states. The output signal y is generated based on C matrix (iutput matrix) multiplied by the states and the D matrix (feedthrough matrix) directly passing through the input signals to the output.
The spreadsheet layout is presented in fig. 1. The top part contains the data you can specify as the user (marked in yellow). This data consists of the A, B, C and D matrices elements, the initial conditions (upper left corner), time step (bottom right corner) and the input signal generator.
Fig. 1. Solving 5x5 MIMO LTI state-space system in a spreadsheet.
The dark yellow cells take values of 0 or 1 and are used to enable and disable plotting of the corresponding input and output channels.
You can solve systems of up to 5 states and use 5 input and 5 output channels as well. The inputs are generated as a sum of 5 sinusoidal functions, which you can parametrize in the panel at the botom of the user input area. The parameters for the signle channel are: the amplitude A, the frequency ω [rad/s], the phase φ [rad] and the starting time t0 [s]:
u_i\;=\; \mathcal{1}(t\:-\:t_0) \: \cdot \: A \: \sin \: (\omega \: t \: + \: \varphi)
The input panel makes it possible to generate a range of popular signals:

Unit step

Unit step is easily made by setting the amplitude A to 1, the frequency ω to 0, and the phase φ to π/2. This effectively makes the generated input to be a cosine calculated at 0 - or simply constant 1.
You can also use the fourth row (t0) to specify at which time the signal should be enabled.

Sine wave

Simply insert the desired values of the amplitude, frequency and phase.

Square wave

We can use Fourier expansion to approximate a square wave of period L as:

f(x) \; = \; \frac{4}{\pi} \: \sum_{n=1,3,5,...} \: \frac{1}{n} \: \sin \: \left( \frac{n\pi x}{L} \right)
An example using 5 input channels to generate a square wave approximation (see fig. 2) is located at examples/square_wave.ods.

Fig. 2. Square wave approximated with Fourier expansion.

Examples

Below are a few examples of the spreadsheet usage.

Integrator

The simplest example is an integrator element, described with the following equation:
y(t) \; = \; \int_0^t \: x(\tau) \: d\tau
Its state-space realization is:

\dot x \; = \; u
y \; = \; \frac{1}{T_i} \: x
where u is the input, x is the integrator state and y is the output. Ti is the integrator time constant - the time at which the integrator output reaches the level of 1 when subjected to an unit step input.
The state-space matrices, based on the equation above are:
A \; = \; [\:0\:]
B \; = \; [\:1\:]
C \; = \; [\:\frac{1}{T_i}\:]
D \; = \; [\:0\:]
The spreadsheet for this example is available in the examples/integration.ods file. The output of the integrator with the unit step input is shown in fig. 3.
Fig. 3. Integrating an unit step (Ti = 50 s).
Try changing the input to a sine function. The output should then take form of the sine integrated - or a cosine. Hopefully the numerical errors won’t accumulate too much!

Inertia

Another classical example is an inertial element. It represents an object that takes a while to reach its set-point value (e.g. a temperature in a heater, rotations of the engine etc.).
It is described with the following equations:

\dot{x}\;=\;u\:-\:\frac{1}{T}\:x
y\;=\;\frac{k}{T}\:x
You only need one state to deal with the inertial element. The matrices look like this:
A \; = \; [\:-\frac{1}{T}\:]
B \; = \; [\:1\:]
C \; = \; [\:\frac{k}{T}\:]
D \; = \; [\:0\:]
Fig. 4 shows the inertial element response. The example is located in examples/inertia.ods.
Fig. 4. Inertial element response to an unit step (k = 1, T = 10 s).

Mass-spring-damper system

Another classic. A harmonic damped oscillator is one of the cornerstones of physics. It is defined by the following equations:
m\:\ddot{x}\;=\;-k\:x\:-\:b\:\dot{x}\:+\:u
y\;=\;x
We have to use two state variables (x1 representing the position and x2 representing the velocity) to reduce the order of that differential equation and arrive at the following representation:
\dot{x}_1\;=\;\frac{1}{m}\:x_2
\dot{x}_2\;=\;-k\:x_1\:-\:b\:x_2\:+\:u
y\;=\;x_1
The state representation matrices are:
A\;=\; \left[ \: 0, \; \frac{1}{m}; \; -k, \; -b\:\right]
B\;=\; \left[ \: 0; \; 1\:\right]
C\;=\; \left[ \: 1, \; 0 \:\right]
D\;=\; \left[ \: 0 \:\right]
The system response in shown in fig. 5. The example is to be found in examples/mass.ods file.
Fig. 5. Mass-spring-damper system response to an unit step (parameters: m = 1, k = 1, b = 2.5). The blue line is the mass position and the red one is the velocity.

PID

The ideal PID controller response is described by the following differential equation:
y\:(t) \; = \; k_p \: \left(\: e\:(t) \: + \: \frac{1}{T_i}\:\int e\:(\tau) \: d\tau  \: + \: T_d \: \frac{d}{dt} \: e\:(t) \right)
It is not unfortunately possible to implement an ideal derivative, and thus the following approximation is used (here in the form of a transfer function):
G_{PID}\:(s) \; = \; k_p \: \left( \: 1 \: + \: \frac{1}{T_i \: s} \: + \: \frac{T_d \: s}{T \: s \: + \: 1}\right)
Such a PID controller could be realized in the state-space in the following form:
\dot{x}_1 \; = \; u
\dot{x}_2 \; = \; -\frac{1}{T}\: x_2 \: + \: u
y \; = \; \frac{k_p}{T_i}\: x_1 \: + \: \frac{k_p\:T_d}{T} \: x_2  + (k_p \: + \: \frac{k_p\:T_d}{T}) \: u
so that the state-space matrices are:
A\;=\; \left[ \: 0, \; 0; \; 0, \; 1\:\right]
B\;=\; \left[ \: 1; \; 1\:\right]
C\;=\; \left[ \: \frac{k_p}{T_i}, \; -\frac{k_p\:T_d}{T} \:\right]
D\;=\; \left[ \: k_p+\frac{k_p\:T_d}{T} \:\right]
The x1 variable handles the integrator state, while the x2 contains the derivative’s inertial element state. The proportional gain is implemented with the use of the feedthrough matrix D. A sample PID response simulation is presented in fig. 6. The PID controller example can be found in the examples/pid.ods spreadsheet.
Too bad it is somewhat tricky to use spreadsheets for online control!
Fig. 6. Simulating PID controller response (Kp = 1, Ti = 100, Td = 1, T = 1).
It should be easy enough to use the spreadsheet to simulate the behaviour of other models. One idea could be to use it for the simulation of other models I have presented here yet: the cart system, the tanks or the oscillating pendulum (some of these models would have to be linearized first though).