Physical modeling with Math::Model and Perl 6

Let's talk about physics. I'm a physicist by education and profession, so it might be not be too surprising. But I also want to talk about programming, so please bear with me.

When physicists try to understand a system, typically they first come up with a lot of equations, and then try to solve them. Unfortunately some of those equations are quite hard to solve. I've written a module that solves them numerically, and allows you to express the equations in very simple Perl 6 code.

Getting started

Let's start with a very simple model. As a prerequisite you need the Rakudo Perl 6 compiler (latest release or current development version), and then use proto to proto install Math-Model (which also installs its dependencies).

A first model: throwing a stone horizontally into the air

Let's start with an example (to be found in examples/vertical-throw.pl in the Math-Model repository, with small modifications to unconfuse the syntax hilighter used for this page).

It describes the path of stone thrown vertically into the air.

use v6;
use Math::Model;

my $m = Math::Model.new(
    derivatives => {
        y_velocity      => 'y',
        y_acceleration  => 'y_velocity',
    },
    variables   => {
        y_acceleration  => { $:force / $:mass },
        mass            => { 1 },           # kg
        force           => { -9.81 },       # N = kg m / s**2
    },
    initials    => {
        y               => 0,               # m
        y_velocity      => 20,              # m/s
    },
    captures    => ('y', 'y_velocity'),
);

$m.integrate(:from(0), :to(4.2), :min-resolution(0.2));
$m.render-svg('throw-vertically.svg', :title('vertical throwing'));

First we declare that we use Perl 6 (use v6;), and then we load the module that does all the hard work, Math::Model.

Then comes the interesting part: the actual model. It starts by declaring that the derivative of y (which is the name we use for the current height of the stone) is called y_velocity. Likewise the derivative of y_velocity is called y_acceleration.

A derivative describes the rate of change of some other variables. Here is a short list of useful physical quantities and their derivatives:

Next in the model are formulas for some variables. We don't need to give formulas for those values on the right-hand side of the derivatives declarations (y and y_velocity), which we will call integration variables. We just need formulas for the derivatives that are not also integration variables, and for other variables we use in the formulas.

You might remember from physics class that F = m * a, force is mass times acceleration. We use this formula for calculating the acceleration

        y_acceleration  => { $:force / $:mass },

The other variables are actually constants (but Math::Model doesn't distinguish those); mass is set to 1 kg, and the force to the -9.81 N that a mass of 1 kg experiences at the latitude where I live (Germany).

The programmer in you might ask what the heck the : in $:force means. It's a variable that is automatically a named parameter to the current block. Math::Model uses Signature introspection to find out which variables such a block depends on, and topologically sorts all variables into an order of execution that satisfies all the dependencies.

Back to the physical model; all integration variables need an initial value, which are provided on the next few lines. The line

    captures    => ('y', 'y_velocity'),

tells Math::Model which variables to record while the simulation is running.

$m.integrate(:from(0), :to(4.2), :min-resolution(0.2));
$m.render-svg('throw-vertically.svg', :title('vertical throwing'));

Finally these two lines are responsible for running the actual simulation (for the time t = 0s to t = 4.2s), and then write the resulting curves into the file throw-vertically.svg. And here it is, without further ado:

The blue line is the height (in meter), and the yellow line the velocity (in meter/second). Positive velocities mean that the stone is moving upwards, negative that it's moving downwards. After about 4 seconds the stone has the height with which it started, at a velocity of about -20m/s. You don't see any effect of the stone hitting the ground, because the model knows nothing about a ground located at a height of 0m. Just imagine the stone fell into a well or down a cliff, and thus can have negative height too.

Iterating: throwing a stone at an arbitrary angle

It gets a bit more involved when you consider not only the y direction, but also the x direction, and throw the stone at an arbitrary angle. Still the x and y axis are rather independent:

use v6;
use Math::Model;

for 30, 45, 60 -> $angle {

    my $m = Math::Model.new(
        derivatives => {
            y_velocity      => 'y',
            y_acceleration  => 'y_velocity',
            x_velocity      => 'x',
        },
        variables   => {
            y_acceleration  => { $:force / $:mass },
            mass            => { 1 },           # kg
            force           => { -9.81 },       # N = kg m / s**2
            x_velocity      => { 20 * cos($angle, Degrees) } # m / s
        },
        initials    => {
            y               => 0,               # m
            y_velocity      => 20 * sin($angle, Degrees),    # m/s
            x               => 0,               # m
        },
        captures    => <y x>,
    );

    $m.integrate(:from(0), :to(2 * 20 * sin($angle, Degrees) / 9.81), :min-resolution(0.2));
    $m.render-svg("throw-angle-$angle.svg", :x-axis<x>,
                  :width(300), :height(200),
                  :title("Throwing at $angle degreees"));
}

The x velocity stays approximately constant (we disregard air drag), and it's value is the cosine of the throwing angle times the initial, total velocity. The velocity in y direction is the sine of the throwing angle times the initial, total velocity.

As a trick to get nicer graphs in the end, I calculated the expected time until it hits the ground (which I know to be 2 * v_y / a; this is cheating, but it doesn't change the physics behind it).

Finally this time we don't want to have the elapsed time on the x axis, but the actual position in x direction. That can be done with the :x-axis<x> option.

(Note that I cheated to produces these charts; normally SVG::Plot, the module that produces these graphs, automatically scales the axis; I've manually suppressed that; if you try to out the example you'll see that the three charts look much the same, except for different axis scaling).

The maximum reach in x direction before hitting the ground is achieved at 45 degrees. For smaller angles it's not long enough in the air, and for larger angles too much of the velocity is "wasted" in the vertical direction, so the stone doesn't get very far.

A third model: oscillating spring

Finally I want to present a third, simple model: a mass hanging down from a spring, this time including air drag. It doesn't show off any new features of Math::Model, but the result is a nice, oscillating curve.

use v6;
use Math::Model;

my $m = Math::Model.new(
    derivatives => {
        velocity     => 'height',
        acceleration => 'velocity',
    },
    variables   => {
        acceleration    => { $:force / $:mass },
        mass            => { 1 },
        force           => { - $:height - 0.2 * $:velocity * abs($:velocity)},
    },
    initials    => {
        height      => 1,
        velocity    => 0,
    },
    captures    => <height>,
);

$m.integrate(:from(0), :to(20), :min-resolution(1));
$m.render-svg('spring.svg', :title('Spring with damping'));

The air drag/damping is built into the formula for the force. It is proportional to the square of the velocity (making analytical solutions a nuisance), and always points in the opposite direction as the velocity.

Moving on

If you want to have some fun with Math::Model, here's a nice idea: You can simulate a mass attached to a spring, which can move freely in the xy plane. If you chose random initial values for the velocities in both directions, you can get non-periodic (chaotic?) trajectories in the xy plane.

(If you do, please send me the model so that I can also play with it :-))

A personal note

Ever since I've used the simulation tool "stella" in school (must have been year 2002 or so), I've wanted to write something similar. Math::Model is it. It was quite fun to implement in Perl 6. The only drawback is that it's quite slow.