We apply a large hammer to a small problem to illustrate the hammer.

We have two heat reservoirs 0 and 1, a channel between them of conductance P, and another channel of conductance Q between reservoir 1 and constant temperature 0. (Conductance is the reciprocal of resistance.) The rate of heat flow for the first channel is P times the temperature difference T2 − T1. The heat capacity of the reservoir j is Cj. The total heat in reservoir j is TjCj.

C0dT0/dt = P(T1 − T0)
C1dT1/dt = P(T0 − T1) − QT1
When we get to actually computing we will assume that P = 2, Q = 1, C0=10, C1=3. When t=0 T0 = 8 and T1 = 14.

This falls into the more general form
CidTi/dt = ΣjcijTj, where c00 = −P; c01 = P; c10 = P; c11 = −P−Q.
and the naïve code to compute the answer.

Letting Mij = cij/Ci we can rewrite:
dTi/dt = ΣjMijTj

Rewriting in Matrix form we get

dT/dt = MT
Here is a curious scheme that serves for some purposes but is not what we want to study here.

Having invoked matrix notions, we go a step further and invoke vector spaces. We view our set of temperatures T as a vector whose components are the temperatures of the individual components, Our equation tells us how T changes with time, depending on T itself. If we found another coördinate system in this space in which the matrix M' for our transformation was diagonal then in those coördinates we would have n independent equations each of the form dT'i/dt = DiT'i where T'i and Di are real numbers and Di is the ith element on the diagonal of M'. Finding such a transformation is the classic eigenvector problem and the Di’s are the eigenvalues. The matrix for the new coördinate system is M' = G−1MG where G is a matrix whose columns are the eigenvectors.

Our new coördinate system uses the eigenvectors as a basis; Fj is the jth basis vector of the new coördinate system and we can write any solution to these new differential equations as ΣjajFjeDjt for some appropriate set of aj’s. This Algol 68 program finds G and prints G−1MG: a diagonal matrix. The program id “g” denotes G which names a row of column vectors each of which is an eigenvector. We speak of a new coördinate system whose basis vectors are these eigenvectors. Choose an eigenvector x; expressed in the new coördinates it is <0, … 0, 1, 0, … 0>. Gx names the same vector in the old coördinates. MGx names the vector (λx) stretched in the same direction, expressed in the old coördinates. G−1MGx names the stretched vector in the new coördinates. G−1MG should thus be diagonal! That is indeed the last matrix printed and it is diagonal with the eigenvalues on the diagonal!

After printing G−1MG the values ΣjajFjeDjt are printed in the old coördinates, for successive times. The aj are chosen to fit the initial conditions. The answers agree with the simpler calculations.


This C code captures the business end of the above. It omits values calculated only to verify the math. It relies on the routine slv here and routine eigen here.

The line below main() has all the inputs.
In the equation CidTi/dt = ΣjcijTj, “c” is c of the program. The d’s of the program are the reciprocals of the “C”’s in the equation. The definition of n must also be adjusted. sl[j][n] is the aj in the description of the logic above.


Here is an alternate derivation of some of the logic.