Quick and Dirty Complex Functions of Time with Guile + Gnuplot

f(t) = e^(st) where s = 0+i

Inspired by a great video on the Laplace transform, I was looking for an easy way to visualize complex functions of time, i.e., functions of time where the output is a complex number. This is fairly easy to do with Guile Scheme and Gnuplot, though the method could be polished more into a proper application, which I didn’t do.

This is the core function for plotting the data:

(define (complex-plot t0 tmax step f)
  (unless (> t0 tmax)
    (let ((z (f t0)))
      (display t0)
      (display " ")
      (display (real-part z))
      (display " ")
      (display (imag-part z))
      (display "\r\n")
      (complex-plot (+ t0 step) tmax step f))))

We just pass in the starting value (t0), the maximum time value, a step value (usually some fraction of Pi), and a lambda with our complex function of time. Writing complex functions of time is pretty easy in Guile Scheme because all the numbers are represented internally as complex numbers anyway. E.g.:

scheme@(guile-user)> (complex-plot 0 25 (/ 3.1415 32) (lambda (t) (exp (* 0+i t))))
0 1.0 0.0
0.098171875 0.9951850104692727 0.09801425884672994
0.19634375 0.9807864101254524 0.19508464243304188
0.294515625 0.9569428571883648 0.2902763649975122
0.3926875 0.9238839645735444 0.3826727322450213

We have to output this to a file:

(with-output-to-file "out.txt" (lambda () (complex-plot 0 25 (/ 3.1415 32) (lambda (t) (exp (* 0+i t))))))

Then we start gnuplot and use a 3d splot with lines. (Credit goes to this tutorial.)

gnuplot> set xrange [0:25]
gnuplot> set yrange [-2:2]
gnuplot> set ticslevel 0
gnuplot> splot "out.txt" u 1:2:3 with lines

Another interesting complex function of time is the addition of two complex functions of time which leave only a signal on the real axis:

(lambda (t) (+ (exp (* 0+i t)) (exp (* 0-i t))))
e^(s_0 * t) + e^(s_1 * t) where s_0 = 0+i and s_1 = 0-i

And here is my hacker emacs screen, where the action is happening!

hacker emacs screenshot


Representation of Circles by Hermitian Matrices

I borrowed the book Geometry of Complex Numbers by Schwerdtfeger. The material has been fascinating so far, though admittedly it took me about 3 hours to comprehend the material in page 1 of chapter 1. The first subject is this intriguing idea that you can represent a circle as a matrix. More specifically, you you can represent a circle with (complex) center 𝛾 and radius 𝜌 as a matrix

\begin{bmatrix} A & B \\ C & D \end{bmatrix}

Where B = – A𝛾̅, C = -A𝛾, and D = A(𝛾𝛾̅-𝜌²), so that A and D are real numbers, and B and C are complex numbers. For normal circles, you will have A=1, though you can scale the matrix to have other (non-zero) values of A and still have the same circle. (If A is zero, you end up with a straight line, or some other mysterious thing that is not a circle.)

The simple case is the unit circle:

Which is \begin{bmatrix} 1 & 0 \\ 0 & -1 \end{bmatrix}

More exotic is the circle away from the origin, centered at 10+10i, i.e., (10, 10), with radius 3:

\begin{bmatrix} 1 & -10+10i \\ -10-10i & 191 \end{bmatrix}

Here is some code to generate the matrix in Racket, and to put it back:

lang racket

(require math/array)
(require math/matrix)

(define (circle-to-matrix zC r)
  (let ([B (* -1 (conjugate zC))]
        [mzC (magnitude zC)])
    (matrix [[ 1             B                       ]
             [ (conjugate B) (- (* mzC mzC) (* r r)) ]])))

(define (matrix-circle-radius M)
  (let ([A (array-ref M #[0 0])]
        [d (matrix-determinant M)])
    (sqrt (/ d (* -1 (* A A))))))

(define (matrix-circle-center M)
  (let ([C (array-ref M #[1 0])]
        [A (array-ref M #[0 0])])
    (/ C (* -1 A))))
racket@circle.rkt> (circle-to-matrix 0 1)
(array #[#[1 0] #[0 -1]])
racket@circle.rkt> (matrix-circle-radius (circle-to-matrix 0 1))
racket@circle.rkt> (matrix-circle-center (circle-to-matrix 0 1))

Why is this significant? I think Schwerdtfeger is going to cover that on page 2. :)

Fractal with Barycentric Coordinates

This fractal idea did not originate with me, but I wrote some racket code to do the midpoint calculation using barycentric coordinates. This fractal draws a circle at the midpoint of a triangle, then subdivides the triangle and repeats:

Here is the same fractal to four iterations:

To get the midpoints, I could simple pass in the coordinates of the last triangle ABC, and then use “0.5” barycentric coordinates:

        [P1 (barycentric->complex 0.5 0.5 0 A B C)]
        [P2 (barycentric->complex 0.5 0.0 0.5 A B C)]
        [P3 (barycentric->complex 0 0.5 0.5 A B C)]

Here is the full code: