Mon, 10 Dec 2007
To do that I decided I would need a function to calculate the square root of a power series, which I did figure out; it's in the earlier article. But then I got distracted with other issues, and then folks wrote to me with several ways to solve the differential equation, and I spent a lot of time writing that up, and I didn't get back to the original problem until today, when I had to attend the weekly staff meeting. I get a lot of math work done during that meeting.
At least one person wrote to ask me for the Haskell code for the power series calculations, so here's that first off.
A power series a0 + a1x + a2x2 + a3x3 + ... is represented as a (probably infinite) list of numbers [a0, a1, a2, ...]. If the list is finite, the missing terms are assumed to be all 0.
The following operators perform arithmetic on functions:
-- add functions a and b add  b = b add a  = a add (a:a') (b:b') = (a+b) : add a' b' -- multiply functions a and b mul  _ =  mul _  =  mul (a:a') (b:b') = (a*b) : add (add (scale a b') (scale b a')) (0 : mul a' b') -- termwise multiplication of two series mul2 = zipWith (*) -- multiply constant a by function b scale a b = mul2 (cycle [a]) b neg a = scale (-1) aAnd there are a bunch of other useful utilities:
-- 0, 1, 2, 3, ... iota = 0 : zipWith (+) (cycle ) iota -- 1, 1/2, 1/3, 1/4, ... iotaR = map (1/) (tail iota) -- derivative of function a deriv a = tail (mul2 iota a) -- integral of function a -- c is the constant of integration int c a = c : (mul2 iotaR a) -- square of function f sqr f = mul f f -- constant function con c = c : cycle  one = con 1
The really interesting operators perform division and evolve square roots of functions. I discussed how these work in the earlier article. The reciprocal operation is well-known; it appears in Structure and Interpretation of Computer Programs, Higher-Order Perl, and I presume elsewhere. I haven't seen the square root extractor anywhere else, but I'm sure that's just because I haven't looked.
-- reciprocal of argument function inv (s0:st) = r where r = r0 : scale (negate r0) (mul r st) r0 = 1/s0 -- divide function a by function b squot a b = mul a (inv b) -- square root of argument function srt (s0:s) = r where r = r0 : (squot s (add [r0] r)) r0 = sqrt(s0)We can define the cosine function as follows:
coss = zipWith (*) (cycle [1,0,-1,0]) (map ((1/) . fact) [0..])We could define the sine function analogously, or we can say that sin(x) = √(1 - cos2(x)):
sins = (srt . (add one) . neg . sqr) cossThis works fine.
f = f0 : mul2 iotaR (deriv f)This holds for any function, and so it's unsolvable. But if you combine it with the differential equation, which says that f = f', you get:
f = f0 : mul2 iotaR f where f0 = 1 -- or whatever the initial conditions dictateand in fact this works just fine. And then you can observe that this is just the definition of int; replacing the definition with the name, we have:
f = int f0 f where f0 = 1 -- or whateverThis runs too, and calculates the power series for the exponential function, as it should. It's also transparently obvious, and makes me wonder why it took me so long to find. But I was looking for solutions of the form:
f = deriv fwhich Haskell can't figure out. It's funny that it only handles differential equations when they're expressed as integral equations. I need to meditate on that some more.
It occurs to me just now that the f = f0 : mul2 iotaR (deriv f) identity above just says that the integral and derivative operators are inverses. These things are always so simple in hindsight.
Anyway, moving along, back to the original problem, instead of f = f', I want f2 + (f')2 = 1, or equivalently f' = √(1 - f2). So I take the derivative-integral identity as before:
f = int f0 (deriv f)and put in √(1 - f2) for deriv f:
f = int f0 ((srt . (add one) . neg . sqr) f) where f0 = sqrt 0.5 -- or whateverAnd now I am done; Haskell cheerfully generates the power series expansion for f for any given initial condition. (The parameter f0 is precisely the desired value of f(0).) For example, when f(0) = √(1/2), as above, the calculated terms show the function to be exactly √(1/2)·(sin(x) + cos(x)); when f(0) = 0, the output terms are exactly those of sin(x). When f(0) = 1, the output blows up and comes out as [1, 0, NaN, NaN, ...]. I'm not quite sure why yet, but I suspect it has something to do with there being two different solutions that both have f(0) = 1.
Sample code is here. For a Scheme implementation, see SICP. For a Java, Common Lisp, Python, Ruby, or SML implementation, do the obvious thing.
But anyway, it does work, and I thought it might be nice to blog about something I actually pursued to completion for a change. Also I was afraid that after my week of posts about Perl syntax, differential equations, electromagnetism, Unix kernel internals, and paint chips in the shape of Austria, the readers of Planet Haskell, where my blog has recently been syndicated, were going to storm my house with torches and pitchforks. This article should mollify them for a time, I hope.
[ Addendum 20071211: Some additional notes about this. ]