## Computing the steady state using Tellurium

If you’re building a model and you want to quickly find the model’s steady state, you can call the command steadyState. Let’s illustrate this with an example:

import tellurium as te

# Define a simple linear chain of reactions
$Xo -> S1; k1*Xo; S1 -> S2; k2*S1; S2 -> S3; k3*S2; S3 ->$X1; k4*S3;

# Initialize rate constants
k1 = 0.2; k2 = 0.7; k3 = 0.15; k4 = 1.3;
Xo = 10
''')



Running this script by clicking in the green arrow in the tool bar will output the following steady state levels:

[2.85714286  13.33333333   1.53846154]


But which steady state value refers to which species in my model? To find that out call getFloatingSpeciesIds:

print r.getFloatingSpeciesIds()
['S1', 'S2', 'S3']


The order of the speces Ids match the order of steady state values. In other words, S1 = 2.857, S2 = 13.333, and S3 = 1.538.

If you want to check that these values really are the steady state values you can compute the rates of change:

print r.getRatesOfChange()
array([  0.00000000e+00,  -4.44089210e-16,   4.44089210e-16])


If you look carefully all the rates of change are effectively zero, thus confirming we’re at steady state.

What about the stability of the steady state? That is, is it stable to small perturbations in S1, S2 or S3? To find this out we need to compute the eigenvalues of the system Jacobian matrix. If all the eigenvalues are negative this means small pertubrations are damped so that the system returns to the steady state. There are two ways to do this, we can use the getReducedEigenValues call, or we get compute the Jacobian first and then compute the eigenvalues of the Jacobian. Both these approaches are given below. Its probably simplest to call just getReducedEigenValues.

print r.getReducedEigenValues()
[-0.7  -0.15 -1.3 ]


Notice that all the eigenvalues are negative indicating that small perturbations are stable.

This is the alternative approach which computes the Jacobian first and then the eigenvalues of the Jacobian:

print te.getEigenvalues(r.getFullJacobian())
[-0.7  -0.15 -1.3 ]


Note that we’re using a utility method from the tellurium library, getEigenValues to compute the eigenvalues.

## Bifurcation Analysis with Tellurium

I thought I’d try and write a series of HowTos on Tellurium, our python based tool for the construction, simulation and analysis of biochemical models. Details on this tool can be found here.

One the unique features of Tellurium is that it comes with the AUTO2000 package. This is a well established software library for performing a bifurcation analysis.

Unlike other implementations (not includng oscill8), AUTO2000 in Telluirum does not require a separate compiler to compile the model for AUTO2000 to compute the differential equations. This makes it easier to delpoy and users don’t have to worry about such details. AUTO2000 uses libRoadRunner to access and compute the model.

Let’s begin with an example from my text book “Systems Biology: Introduction to Pathway Modeling”, Figure 12.20, page 279 in revision 1.1 of the book. The model in question is a modified model from the ‘Mutual activation’ model in the review byTyson (Current Opinion in Cell Biology 15:221–231, Figure 1e). In this example increasing the signal results in the system switching to the high state at around 2.0. If we reduce the signal from a high level, we traverse a different steady state. If we assume the signal can never be negative, we will remain at the high steady state even if the signal is reduced to zero. The bifurcation plot in the negative quadrant of the graph is physically inaccessible. This means it is not possible to transition to the low steady state by decreasing signal. As a result, the bistable system is irreversible, that is, once it is switched on, it will always remain on. To compute the bifurcation diagram we first define the model:

import tellurium as te
import teplugins
import pylab as plt

$X -> R1; k1*EP + k2*Signal; R1 ->$w;  k3*R1;
EP -> E;   Vm1*EP/(Km + EP);
E -> EP;   ((Vm2+R1)*E)/(Km + E);
Vm1 = 12; Vm2 = 6;
Km = 0.6;
k1 = 1.6; k2 = 4;
E = 5; EP = 15;
k3 = 3; Signal = 0.1;
''')


We’ve imported three packages, tellurium to load the model, teplugins to access AUTO2000 and pylab to gain acess to matplotlib. Once we have the model loaded we can get a handle on AUTO2000 by calling teplugins.Plugin(“tel_auto2000”) and set a number of properties in AUTO2000. This includes loading the model into AUTO200, identifying the parameter we wish to modify for the bifurcation diagram (in this case signal), following by some options to carry out a pre simulation to help with the initial location of the steady state and finally the limits for x axis for the plot, in this case -2 to 3. Details of other propoerties to change can be found by typing auto.viewManual(), make sure you have a pdf reader available. The alternative is to go to the intro page.

auto = teplugins.Plugin("tel_auto2000")

auto.setProperty("SBML", r.getCurrentSBML())
auto.setProperty("ScanDirection", "Negative")
auto.setProperty("PrincipalContinuationParameter", "Signal")
auto.setProperty("PreSimulation", "True")
auto.setProperty("PreSimulationDuration", 1.0)
auto.setProperty("RL0", -2.0)
auto.setProperty("RL1", 3.0)


To run the bifurcation analysis we use the Python code:

auto.execute()


If all was successful we can next plot the results. It is possible to plot the results using your own code (see below) but it might be more convenient to use the builtin facilties, for example:

pts     = auto.BifurcationPoints
lbls    = auto.BifurcationLabels
biData  = auto.BifurcationData
biData.plotBifurcationDiagram(pts, lbls)


The pts vector contains the point coordinates where the bifurcation points are located. lbls give the labels that correspond to the pts vector and indicate what type of bifurcation point it represented. Finally a special object, here called biData contains the data together with a number of useful utilties. The import important of these is biData.plotBifurcationDiagram(pts, lbls) which takes pts and lbls as arguments. Running this code will generate the bifurcation plot shown below.

We can also print out a text summary of the computation using the command, auto.BifurcationSummary, which returns a summary of tjhe fidnnngs.

Summary :
BR PT TY LAB PAR(0) L2-NORM U(1) U(2)
1 1 EP 1 3.00000E+000 2.38908E+001 1.42336E+001 1.91879E+001
1 50 2 5.48632E-001 2.14502E+001 1.06592E+001 1.86144E+001
1 95 LP 3 -9.31897E-001 1.79088E+001 7.44448E+000 1.62881E+001
1 100 4 -9.13592E-001 1.74094E+001 7.22867E+000 1.58377E+001
1 150 5 1.74114E-001 1.26908E+001 6.15211E+000 1.10999E+001
1 200 6 1.54579E+000 8.36095E+000 5.44499E+000 6.34488E+000
1 234 LP 7 2.06965E+000 5.51178E+000 4.47540E+000 3.21723E+000
1 250 8 1.84381E+000 4.03137E+000 3.51306E+000 1.97747E+000
1 300 9 -5.77433E-001 7.02653E-001 -5.14930E-001 4.78089E-001
1 325 EP 10 -2.00625E+000 2.56175E+000 -2.55121E+000 2.32100E-001


We can manually plot the data by gaining access to the numpy version of the data. To do this we use:

pltData = biData.toNumpy


pltData is a numpy array where the first column is the bifurcation parameter and the remaning columns contain the species. For example to plot the bifurcation diagram for the first species in the model, R1 we would use:

plt.plot(x[:,0], x[:,1], linewidth=2)
plt.axvline(0, color='black')
plt.xlabel ("Signal")
plt.ylabel ("R1")


I added a axvline command to draw a vertical line from the zero axis. I also added some axis labeling statements. These commands will result in:

What is interesting about this model is that the upper branch reaches the zero parameter value before the turning point. This means it is difficult to switch to the lower steady state by just lowering the signal.

Viewing thee Network

One other things we can do is view the model as a network. Telluirum comes witha simple network viewer in the package nwed. import the viewer using

import nwed


at the ipython console. To view the network make sure the network viewer panel is visible, do this by going to the View menu, find panes and select, then look down the menu items and near the bottom you’ll find Network Viewer, select this option. To view the network, type the following at trhe ipython console.

nwed.setsbml (r.getSBML())


The viewer should now display something like:

Note that every view will be different and depends on the layout algorithm.

## Tikz Code for Drawing Metabolic Feedback Loops

I needed some figures that displayed a variety of different negative feedback loops so I created these using Tikz. Nothing particularly special. There are some absolute distances in the code which perhaps could be removed to make it more generic.

\documentclass{article}

\usepackage{amsmath}
\usepackage{tikz}
\usetikzlibrary{arrows}
\usetikzlibrary{calc}

\begin{document}

\begin{tikzpicture}[>=latex', node distance=2cm]
\node (Xo) {};
\node [right of = Xo] (S1) {\Large $x$};
\node [right of = S1] (S2) {};

\draw [->,ultra thick,blue] (Xo) -- node[above, black] {$v_1$} (S1);
\draw [->,ultra thick,blue] (S1) -- node[above, black] {$v_2$} (S2);

% Lets draw a line with a blunt end, -|
\draw [-|,ultra thick,blue]
% start in the middle of S2, and move down 2.75 mm
% ($...$) notation is used to add the coordinates
($(S1) + (0mm,-2.75mm)$)

% Now draw the line down by 3mm
% -- means draw to, + means move by
-- +(0,-3mm)
% Now move back to the left of S2
% The symbol -| means draw horizontal then vertical.
% If we used -- instead the line would be drawn
% diagonally to the reaction edge.
% (S1) is the center of the node. But we want the
% blunt end to end below the S1 line and
% half way to the left. The 10mm is half the node
% distance of 2cm, and 1mm is slightly below
% the reaction line.
-| ($(S1) - (10mm,1mm)$);
\end{tikzpicture}

\vspace{1cm}

\begin{tikzpicture}[>=latex', node distance=2cm]
\node (Xo) {};
\node [right of = Xo] (x1) {\Large $x_1$};
\node [right of = x1] (x2) {\Large $x_2$};
\node [right of = x2] (x3) {};

\draw [->,ultra thick,blue] (Xo) -- node[above, black] {$v_1$} (x1);
\draw [->,ultra thick,blue] (x1) -- node[above, black] {$v_2$} (x2);
\draw [->,ultra thick,blue] (x2) -- node[above, black] {$v_3$} (x3);

% Lets draw a line with a blunt end, -|, using the following coords
\draw [-|,ultra thick,blue]
% start in the middle of x2, and move down 2.75 mm
% ($...$) notation is used to add the coordinates
($(x2) + (0mm,-2.75mm)$)

% Now draw the line down by an additional 3mm
% -- means draw to, + means move by
-- +(0,-3mm)

% Now move back to the left of x2
% The symbol -| means draw horizontal then vertical.
% If we used -- instead the line would be drawn
% diagonally to the reaction edge.
% (S1) is the center of the node. But we want the
% blunt end to end below the S1 line and
% half way to the left. The 10mm is half the node
% distance of 2cm, and 1mm is slightly below
% the reaction line.
-| ($(x1) - (10mm,1mm)$);
\end{tikzpicture}

\vspace{1cm}

\begin{tikzpicture}[>=latex', node distance=2cm]
\node (Xo) {};
\node [right of = Xo] (x1) {\Large $x_1$};
\node [right of = x1] (x2) {\Large $x_2$};
\node [right of = x2] (x3) {\Large $x_3$};
\node [right of = x3] (x4) {};

\draw [->,ultra thick,blue] (Xo) -- node[above, black] {$v_1$} (x1);
\draw [->,ultra thick,blue] (x1) -- node[above, black] {$v_2$} (x2);
\draw [->,ultra thick,blue] (x2) -- node[above, black] {$v_3$} (x3);
\draw [->,ultra thick,blue] (x3) -- node[above, black] {$v_4$} (x4);

% Lets draw a line with a blunt end, -|
\draw [-|,ultra thick,blue]
($(x3) + (0mm,-2.75mm)$)
-- +(0,-3mm)
-| ($(x1) - (10mm,1mm)$);
\end{tikzpicture}

\vspace{1cm}

\begin{tikzpicture}[>=latex', node distance=2cm]
\node (Xo) {};
\node [right of = Xo] (x1) {\Large $x_1$};
\node [right of = x1] (x2) {\Large $x_2$};
\node [right of = x2] (x3) {\Large $x_3$};
\node [right of = x3] (x4) {\Large $x_4$};
\node [right of = x4] (x5) {};

\draw [->,ultra thick,blue] (Xo) -- node[above, black] {$v_1$} (x1);
\draw [->,ultra thick,blue] (x1) -- node[above, black] {$v_2$} (x2);
\draw [->,ultra thick,blue] (x2) -- node[above, black] {$v_3$} (x3);
\draw [->,ultra thick,blue] (x3) -- node[above, black] {$v_4$} (x4);
\draw [->,ultra thick,blue] (x4) -- node[above, black] {$v_5$} (x5);

% Lets draw a line with a blunt end, -|
\draw [-|,ultra thick,blue]
($(x4) + (0mm,-2.75mm)$)
-- +(0,-3mm)
-| ($(x1) - (10mm,1mm)$);
\end{tikzpicture}

\vspace{1cm}

\begin{tikzpicture}[>=latex', node distance=2cm]
\node (Xo) {};
\node [right of = Xo] (x1) {\Large $x_1$};
\node [right of = x1] (x2) {\Large $x_2$};
\node [right of = x2] (x3) {\Large $x_3$};
\node [right of = x3] (x4) {};

\draw [->,ultra thick,blue] (Xo) -- node[above, black] {$v_1$} (x1);
\draw [->,ultra thick,blue] (x1) -- node[above, black] {$v_2$} (x2);
\draw [->,ultra thick,blue] (x2) -- node[above, black] {$v_3$} (x3);
\draw [->,ultra thick,blue] (x3) -- node[above, black] {$v_4$} (x4);

% Lets draw a line with a blunt end, -|
\draw [-|,ultra thick,blue] (10mm, -8mm)
-- +(0,6mm);
\node (x) at (10mm,-11mm) {\Large $x$};
\end{tikzpicture}

\end{document}


## The Confusion of Modern Textbooks – Stylistic Sugar

I’ve been looking for a textbook on statistics for a class I’ll teach in the autumn term. While the content of many textbooks might be ok the way the information is presented makes them difficult to read – at least it does for me. This applies to most undergraduate textbooks published today whether they be about statistics or other topics. There seems to be a need by publishers to embellish textbooks with so much stylistic sugar that the content is buried.

I scanned two typical pages from a second-hand stats textbook I bought from Goodwill to illustrate what I mean.

I’ve marked using a red star the many different styles used in the text book on two random pages, these include:

1. Section heading with a vertical dark purple line
2. Highlighted area using three colors (black, blue and purple)
3. A notes section with blue heading
4. An Illustration section using blue font but spaced out letters and thin left-bar in dark purple
4. Up and down purple arrows in illustration section indicating question and answer
6. Paragraph of text using back font (finally something normal)
7. Figure caption in dark purple text with green vertical line and horizonal line in black
8. Exercise box in margin with heading in white with black background
9. Exercise box where the question in black with yellow background
10. Typewriter font for computer code with black and blue horizontal lines delimiting code

Not included in these pages are also other stylistic sugar:

1. Case study section that uses five colors and six stylistic features
2. Exercises with six stylistic features including at least eleven different symbols
3. Call-out in exercises using a script font with blue background.
3. Chapter practice tests in black font with blue background and thick blue horizontal line
4. Chapter Objective in black font, red bullet point in off yellow background with vertical dotted line.
5. Chapter heading page, eight stylistic features with multiple fonts

And this is before a student has even started to read the content I counted at least 15 different fonts used in the text.

## Transistor Based Flip-Flop or a (not)RS NAND Latch

Been a while since I did a posting, too much time spent writing grants, and I mean a lot of time. In this blog I thought I’d describe a small project I did a month or two ago to build a sinple RS NAND latch using transitors only.The RS stands for Reset/Set. The latch itself is based on two NAND gates connected to each other in a cycle. The following diagram (borrowed from Wikipedia) shows the latch in terms of two NAND gates:

The inputs are designated bar S and bar R (bar meaning not). The outputs are Q and bar Q. Let us assume bar S is set to digital input zero and bar R to one. If you follow the logic through the two NAND gates you’ll realise that this means that Q will have a value of zero and bar Q one. This is the reset state and relative to Q, the circut stores a value of zero. If we now set bar S to one and bar R to zero, the circuit will flip to a new state where Q is now at one and bar Q at zero. Relative to Q the circuit stores a value of one. The point of the circuit is that even if bar S (or bar R) now goes back to zero the circuit remains in its last state. The only way to change the circuit is to apply a one to the reset or set input lines. The two states that the circuit exhibits are stable.

One can buy flip-flops ready made and one can also buy NAND gates ready made (7400). What I wanted to do was build a flip-flop using discrete transistors. I found the following site that describes how to build a 4 bit adder using transitors and from their circuit diagrams I came up with the following NAND circuit using transitors:

In the orgiinal circuit the input resistors were 10K but I wasn’t able to get enough forward bias with 10K so I dropped the input resistors to 4.7K. As with the adder circuit I used common NPN BC 547 tranistors. They are pretty cheap on ebay, you can get 50 for a dollar (!). I used 10 volts for the voltage supply but 9 volts willl work too. I used the following PCB boards, cut in half, to make a single NAND circuit:

Again you can find these cheaply on ebay for about 7 dollars for 10 boards. The image below shows a single NAND gate made using the circuit above.

I made two of the above NAND gates andf connected them together. I also made two LED driver circuits that would cature the outputs Q and bar Q. The final circuit looks like:

There are two push buttons on the broad board next to the Set and Reset labels which momentarily take the inputs high (one), turning on the transitors. The circuit on the right that uses the white PC board is a simple LED driver, this also uses a BC547, with a 1K resistor on the output and a 4.7K resistor on the input. When the output from the flipflop goes to one current flows into the LED driver transistor base which in turns switches on the transitor allowing current to flow through LED which lights up. The LED circuit is shown below:

The following video shows the circuit in operation. For those who might be interested in taking this idea to the extreme, I suggest you check out the megaprocessor.

## Z80 Microcomputer

Its been a while since I wrote something for this blog, mainly due to pressures of work. However its summer now and I’ve managed to carve out some time to do other work related projects. One thing I’ve done in one of my undergraduate classes is host an honors projects which students can opt into. Part of this is to prove to me they can soldering electronic components onto a PCB board. Every engineer should be able to build electronic gadgets from scratch and solder is an essential skill to have.

In order to keep myself in hardware practice this summer I decided to go back in time and search out a Z80 microprocessor kit to build. When I was a teenager I built and designed my own Z80 system. If I remember right it had 512 bytes of RAM, switches and LED lights to enter and run programs. Once the machine was switched off however it forget any program that had been entered. I’d planned to add a hex keyboard etc but my money and spare time ran out. If you search the web you’ll find many Z80 based projects in various degrees of sophistication. There is an entire community of Z80 builders. Some specific ones that caught me eye include:

Those too young to know what a Z80 is, its one of the early microprocessors that came out in the 1970s, 1976 to be specific. It was basically an enhanced version of the 8080 a very early and popular Intel microprocessor. What was key about the Z80 was that is was quite easy to construct a basic computer using the Z80 chip. The chip required very few support chips, for example the clock was a simple single-phase clock (easily made from two NAND or NOT gates) and it only required a single 5v power supply. These and other features made the Z80 a very popular microprocessor at that time. Many famous and not so famous pre-IBM PC computers were build using the Z80, including the TRS-80, Nascom 1 and 2 (I had the Nascom 2), the famous Sinclair ZX80, ZX80 and ZX Spectrum, Amstrad CPC, Jupiter Ace, and many others. You can still purchase the chip new for only a few dollars.

My problem is I don’t have the time to design a new system for myself so instead, I hunted out any available kits. I found one and a nice one at that called the CPUville Z80 computer kit. This kit was developed and is sold by Donn Stewart. He has various addons to the base CPU board such as a serial interface and disk and memory expansion. What is also nice is that he adds the facility to single clock the microprocessor at very low speeds so that one can observe the operation of the computer in detail. Given my birthday is in July I decide to order one of his kits for myself. His kits are very modestly priced, the basic CPU board (I fully working Z80 computer) only costs \$42 dollars. The kit comes fully complete except for the 5v power supply, but Donn can also supply this. The computer has 2K of RAM and a preprogrammed 2K ROM with various test routines included. Data is entered through a bank of 16 switches and output is via a set of 16 LEDs. To program this computer you’ll need to know binary and hex, and no you can’t browse the net using this computer.

Having obtained the kit, I assembled the various components. The kits comes with an excellent assembly and operation manual. The kits itself only supplies a few chip sockets for the CPU, RAM and ROM. I tend to be more cautious so I purchased a bunch of 20, 8 and 7 pin chip sockets. It probably took me about 6 hours to build. The board has about 500 soldered connections so you have to take your time and check the quality of your work. I took some pictures of the board during construction and a video to show the computer working. Amazingly it worked first time. Operation involves flipping a switch to halt the CPU, entering an address in ROM where a test program is lcoated and then releasing the CPU to run. The test program I ran was a simple counter than counted from zero to 65536 and output the binary number to the LED lights. If you need some fun with a very basic computer this seems to be the one to get.

## Bode Plots using Python

I needed a quick way to plot some Bode plots for a second order system. I didn’t have access to Matlab, instead I searched for a solution using Python, and I found one. Documentation is a bit sparse  so this example might be helpful.

The signals packages supports the signal.bode method which turned out to be quite easy to use. Signal is part of the scipy package and is something we bundle with our Tellurium platform.

There a more comprehensive discussion of Python and Control Theory here.

#
import numpy as np
from scipy import signal
from matplotlib import pyplot as plt

# Coefficients in numerator of transfer function
num = [1]
# Coefficients in denominator of transfer function
# High order to low order, eg 1*s^2 + 0.1*s + 1
den = [1, 0.1, 1]

# Scan over zeta, a parameter for a second-order system
zetaRange = np.arange(0.1,1.1,0.1)

f1 = plt.figure()
for i in range(0,9):
den = [1, 2*zetaRange[i], 1]
print den
s1 = signal.lti(num, den)
% Specify our own frequency range: np.arange(0.1, 5, 0.01)
w, mag, phase = signal.bode(s1, np.arange(0.1, 5, 0.01).tolist())
plt.semilogx (w, mag, color="blue", linewidth="1")
plt.xlabel ("Frequency")
plt.ylabel ("Magnitude")
plt.savefig ("c:\\mag.png", dpi=300, format="png")

plt.figure()

for i in range(0,9):
den = [1, zetaRange[i], 1]
s1 = signal.lti(num, den)
w, mag, phase = signal.bode(s1, np.arange(0.1, 10, 0.02).tolist())
plt.semilogx (w, phase, color="red", linewidth="1.1")
plt.xlabel ("Frequency")
plt.ylabel ("Amplitude")
plt.savefig ("c:\\phase.png", dpi=300, format="png")
%


Output from Python script:

Magnitude plot:

Phase plot:

| 1 Comment

## Disappointed with Drobo 5N

I got myself a Drobo 5N to serve as a backup of my data and other work. For those who don’t know the Drobo 5N is a 5 bay network attached storage (NAS) device. What is interesting about the Drobo is that any of the drives (and more than one) in one of the 5 bays can be removed without any effect to your stored data. This is because all data is redundantly spread across all drives. Not only that but the drives you can have in the 5 bays can be any size and any make.

The problem however is that the Drobo 5N is not truly a network attached storage device because you can only access directly from your computer on your local network, ie your office. So for example if you have a Drobo drive in your office you can’t access your data from home or anywhere else in the world. Most NAS system do allow access, for example by running an ftp server. It is possible to install an ftp server on a Drobo but the installation is buggy and trashes the firmware. So all in all the Drobo 5N is a bit of a disappointment.

## Programming Language Popularity

As programmers we sometimes like to play the game of what is the most popular programming language? There are various metrics online that try to measure this, most notably TIOBEPYPL and Trendy Skills. They are all flawed in someway and often diverge considerably in their rankings. For example PYPL ranks Swift at number 9 whereas its not even in the top  24 on TIOBE. C is ranked 8th on Trendy Skills but 1st in TIOBE and 6th on PYPL. Sometimes they correspond for example, Javascript is ranked 7th place on both TIOBE and PYPL,

Another ranking system, also flawed, is the #code20XX twitter based rank. This is where people vote for what programming languages they used in the specified year. On #code2104 this year, Javascript was the winner closely followed by Python, Java then Ruby. On TIOBE Ruby is ranked 15th, and Python 8th. The other language to show itself on code2014 was Object Pascal. I get a lot of criticism for programming in Object Pascal/Delphi but I like it because its easy to read, very flexible to use and through the various GUI frameworks is ideal for writing cross platform GUI applications. On #code2014 Object Pascal/Delphi came 9th (compared to 20th on TIOBE). Object Pascal/Delphi beat a couple of mainstream languages include C and C++. I think this says more about the enthusiasm of the community than anything else which is of course an important factor when considering a language to use. The downside of Object/Pascal/Delphi is the dollar expense of buying into the platform. Embarcardero who sell the compiler and IDE charge a significant amount to be a member of this group. It effectively excludes most students and hobbyists from the community which is a great shame. On the plus side we have of course freepascal and lazarus which are also excellent cross platform tools and largely Delphi compatible. They are still perhaps a little immature in places but they are developing quickly.

## Simple RK4 Code

Thought I’d start off the New Year with a simple freebie, a 4th Order Runge-Kutta class for Delphi, Windows and Mac OS. Should also work on Android and iOS mobile platforms and with Free Pascal. The code below shows how you’d use it. First declare the set of differential equations that need to be solved (TDoubleDynArray type can be found in Types):

 TMySystem = class (TObject)
class procedure func (time : double; y, p, dydt : TDoubleDynArray);
end;

class procedure TMySystem.func (time : double; y, p, dydt : TDoubleDynArray);
begin
dydt[0] := p[0] - p[1]*y[0];
dydt[1] := p[1]*y[0] - p[2]*y[1];
end;


Next, create a Runge-Kutta object:

// First argument = number of ODEs
// Second argument = number of parameters, if any
// Third argument = ODE Function
rk4 := TRK4.Create (2, 3, TMySystem.func);
rk4.p[0] := vo; rk4.p[1] := k1; rk4.p[2] := k2;

The above code include the option to declare and initialize some parameters that are part of the differential equations. To actually integrate the system by one step use the line:

// Return new time point, assign to startTime ready for next time
startTime := rk4.eval (startTime, y, stepSize);


The y argument will contain the updated solution. Run the line again to do the next step etc. Note that the eval method takes three arguments, that current value for the independent variable (time), an array that contains the current values for the dependent variable, and the step size to use.