## Grapher: a hidden treasure

From Mac OS X 10.4 onward, there is a software tool pre-installed called Grapher, hidden deep in /Applications/Utilities. Grapher is a graphing calculator, which allows users to plot many types of 2-D and 3-D graphs from data points, user-defined functions, and even solutions of differential equations (both implicit and explicit, and up to second-order). It can also draw vector fields, very useful for studying differential equations and dynamical systems. However, as a graphing calculator, it only supports plotting graphs of functions or series, not graphic charts (as you will see in Excel).

I call Grapher a hidden treasure for two reasons. First, it is hidden deep in the Mac OS, so not many Mac users know about it, including math teachers and researchers who often need to draw graphs. Second, it is really a treasure for it is very powerful, interactive, and easy to use. If your graphing task is in its capability, it is always easier to use Grapher than other more sophisticated tools, e.g. Matlab, matplotlib, Sagemath. Furthermore, the outcomes are almost always prettier and more satisfactory. In this post, I do not intend to teach how to use Grapher, not even a short tutorial. Rather, I will present a real case in which I used Grapher for my research and presentation, which not only saved me a lot of time but also gave me accurate and beautiful results.

## Task: plot solution of periodic switched system

In my recent research and presentation, I needed to study the state trajectory of a switched 2-dimensional dynamical systems under periodic switchings. I also needed to study its behavior (i.e. how the trajectory would change) when the sampling period changed. The system under study was more complex, but I will use a simple system for illustration. The system is given by the following switched linear differential equation:

$\displaystyle \dot{x}(t) = \begin{cases} \begin{bmatrix} -0.2 & 0 \\ 0.1 & -0.3 \end{bmatrix} x(t) & \text{if } \sigma(t) = 1 \\ \begin{bmatrix} -0.1 & -0.05 \\ 0 & 0.1 \end{bmatrix} x(t) & \text{if } \sigma(t) = 2 \end{cases}$

where $\sigma(\cdot)$ is the switching signal, and the initial condition is $x(0) = [2, 2]^T$. As you can see, there are two discrete modes and in each mode, the dynamics is linear. The second mode is unstable (the second eigenvalue is 0.1).

We want to study the behavior of such system under a periodic switching signal of the form:

$\displaystyle \sigma(t) = \begin{cases} 1 & \text{if } t \text{ mod } T \in [a, b] \\ 2 & \text{otherwise} \end{cases}$

in which $T > 0$ is the sampling period and $0 \leq a < b \leq 1$ define the segment of the period during which the system is in mode 1.

## The code

The following figure displays the code in Grapher. It is very intuitive and follows almost exactly the mathematical descriptions above. Observe that it is easy to define differential equations and piecewise functions in Grapher.

As a comparison, I implemented the same simulation in Matlab. If this were implemented in Python or Sagemath, the code would have been similar. Because Matlab and its ODE solvers do not naturally support piecewise functions, I had to implement the switchings explicitly by stopping the numerical integration exactly at the switching instants then restarting it. There are two approaches that I can think of:

1. You calculate the switching instants explicitly as $t_0=0, t_1=a T, t_2=(a+b)T, \dots$. Then you call the ODE solver from $t_0$ to $t_1$ with $x(t_0) = x(0)$, save the outputs, call it again from $t_1$ to $t_2$ with $x(t_1)$ from the previous result, concatenate the outputs, and repeat. The code is quite cumbersome and unintuitive. I did not implement this approach but the next.
2. You instruct the ODE solver to detect for events which are the switching instants by defining an event function. This way you can implement the differential equation as a piecewise function using a simple if-else statement.

The Matlab code I wrote follows:

function periodic_switched_ode(x0, tfinal, T, a, b)
[tout, yout, te, ie, ye] = ode45(@func, [0, tfinal], x0, odeset('Events', @events));
plot(yout(:,1), yout(:,2));
disp(te);

function dy = func(t,y)
tt = mod(t, T)/T;
if tt >= a && tt < b
dy = [-0.2, 0; 0.1, -0.3] * y;
else
dy = [-0.1, -0.05; 0, 0.1] * y;
end
end

function [value,isterminal,direction] = events(t,y)
tt = mod(t, T)/T;
value = tt - [a; b];
isterminal = [0; 0];
direction = [0; 0];
end
end

As you can see, it sort of resembles the mathematical description of the problem, but not as close and intuitive as the Grapher code.

Update on October 31, 2011: I wrote a Matlab function called ttode to simulate this kind of systems using approach 1. For more information and to download, please visit my GitHub repository.

## Results

Grapher immediately solved the system and returned the following graph:

The quality of the web image is not high, but Grapher allows exporting to EPS and PDF, which are suitable for use in scientific publications.

The same simulation in Matlab gave the following figure:

Oh wait! The graph was totally different from that of Grapher. Our intuition tells us that the result returned by Matlab was wrong. But what caused the error? The printed event time instants were completely off. But why? Using a supposedly more accurate ODE solver like ode113 gave even worse results. The underlying problem was that Matlab used a coarse estimation of time steps for integration, which skipped many events. So in order to force Matlab to use finer time steps, we need to reduce the MaxStep parameter of the solver as following:

[tout, yout, te, ie, ye] = ode45(@func, [0, tfinal], x0, odeset('Events', @events,'MaxStep',T/50));

As you can see, I set the maximum time step to one-fiftieth of the period. If it had been set to T or even T/2, the results would have been still wrong. The general rule is that the maximum time step should be smaller than the smallest interval between two consecutive events. Sound complicated? Yes, because it is indeed complicated. Especially when you compare it to the Grapher code.

For those who are curious, the switched system is stable. If we set $b = 0.2$, it becomes unstable.

Grapher is a wonderful piece of software. It does an extremely excellent job at what it can do. Of course, Matlab is a larger and more powerful generic tool (and so are Sagemath, Mathematica, etc.). However, Grapher is much more intuitive and easier to use. And prettier too. And it just works (oh Apple). You don’t have to mess up with settings and parameters to make it work correctly. If I hadn’t known the solutions a priori, I might have thought the first result returned by Matlab were correct. Which would have been a disaster in research. Oh, did I mention that Grapher can export the mathematical equations to LaTeX code and images (to be included in slides, for example)?

There are shortcomings of Grapher. Naturally, as a graphing tool, it only supports 2-D and 3-D equations. However, it was not designed as a generic scientific computation software in the first place. Grapher does not allow exporting solutions to data files (e.g. CVS files), which limits its compatibility with other scientific tools. Moreover, it is a Mac-only software, so its files cannot be opened in Windows or Linux.

Regardless, I think Grapher is a wonderful tool for people working with mathematics and graphs. A true treasure somehow hidden in the Mac OS.

Grapher is much more powerful than what I have illustrated. If you want to learn more about Grapher, read this guide and this series of blog posts.