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 },
force => { -9.81 },
},
initials => {
y => 0,
y_velocity => 20,
},
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:
| Quantity | Derivative |
|---|
| position | velocity |
| velocity | acceleration |
| acceleration | jerk |
| momentum | force |
| energy | power |
| charge | current |
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 },
force => { -9.81 },
x_velocity => { 20 * cos($angle, Degrees) }
},
initials => {
y => 0,
y_velocity => 20 * sin($angle, Degrees),
x => 0,
},
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.
