Computational algebra

I’ve been studying electromagnetism recently, and consequently been doing more maths-by-hand. Every time I do this, I think about using computer algebra systems to check my working and help me with calculus. But the only computer algebra system I’ve used was back at uni (Maple, I think) – and whilst Mathematica looks whizzy it’s also quite pricey. So I thought I’d give a few ‘free’ ones a go.

I have a particular ‘test’ problem in mind, because I’ve just worked through it by hand whilst reading Feynman Vol 2. It’s calculating the Laplacian of a radially symmetric potential – which involves a nice mixture of partial and total derivatives, chain rules and product rules. Turns out, I do actually remember how to do all this stuff, but having done it by hand it makes a nice concrete example for doing it on computer.

Another motivation is that I hate the imprecise notation that physicists and some mathematician use. What is the derivative of phi(r) with respect to x? My brain says that phi is a unary function with parameter r, so the only thing that can cause it to change is changes in it’s input r. In physics you have to ‘understand’ that this means phi is really a function over R^3 (space) that’s defined like phi(x,y,z) = g(sqrt(x^2+y^2+z^2)) with the helper g(r)=… telling you how things change with distance. And quite often, phi will also quietly become a function of time too. In SICM they fix this problem by using a scheme-based computer algebra system. Here I’m trying to do the same thing but using a more mainstream maths package.

I just ran sagemath via docker, rather than worry about ‘installing it’:

docker run -it sagemath/sagemath:latest

And so we start to tell it about our world …

%display latex

# We'll use g(r) as our "how things change with distance" function, and we'll leave it abstrac

g = function("g")

# Now we can define phi to be a concrete function over all space that delegates to g

phi(x,y,z) = g(sqrt(x^2+y^2+z^2))
phi.differentiate(x,2)
 
# result
# x^2*D[0, 0](g)(sqrt(x^2 + y^2 + z^2))/(x^2 + y^2 + z^2) - x^2*D[0](g)(sqrt(x^2 + y^2 + z^2))/(x^2 + y^2 + z^2)^(3/2) + D[0](g)(sqrt(x^2 + y^2 + z^2))/sqrt(x^2 + y^2 + z^2)

Assuming that D[0] means ‘derivative with respect to 0th arg’ then this is right but boy is it ugly!

Specifically, we’re not giving a short name ‘r’ to the expression sqrt(x^2+y^2+z^2). We can try to cajoule sagemath along:

r = var('r')
phi.differentiate(x,2).subs({sqrt(x^2+y^2+z^2): r})

Sadly that only improves the g(r) expressions, but not the usages in the denominators.

Perhaps we should tell it about r^2 rather than r?

phi.differentiate(x,2).subs({(x^2+y^2+z^2): r^2})

This does affect the denominator, but weirdly leaves lots of sqrt(x^2) terms. Odd.

We can blunder on ahead and complete the Laplacian:

(phi.differentiate(x,2) + phi.differentiate(y,2) + phi.differentiate(z,2)).simplify_full()

.. which yields the right answer, but again in a not-very-simplified form. Still, this is a pretty impressive feat. It takes a lots of pencil + paper to do these three second derivatives, chain rules and everything else. Doing it in one line shows that sagemath is a pretty useful maths assistant!

So far, we’ve just modelled phi as a three-arg function – but sagemath actually knows what scalar and vector fields are, and can already do operations on them!

from sage.manifolds.operators import *
E.<x,y,z> = EuclideanSpace()
g = function("g")
phi = E.scalar_field(g(sqrt(x^2+y^2+z^2)), name='phi')
phi.laplacian().display()

… which gives the same answer as above, but now we’ve got “batteries included” and can calculate div/grad/curl and all that.

So this is now good enough to solve electrostatic problems, where the laplacian of the scalar potential is proportional to charge density. For dynamic cases however, we’ll need a function of x,y,z and also t. Not sure how that’ll work with the manifold support. Do I need a function from time to scalar fields? A field of time-to-scalar functions? A 4d spacetime structure?

This page is working with E/B fields using 4d manifolds. I’ve done courses on special relativity, so I can handle spacetime but the language of differential geometry (manifolds, 2-forms) is beyond me currently (although Sussman of SICP fame has a similarly computational approach to it). Perhaps it’s better to leave the sage.manifolds alone for now and return to explicit x/y/z/t approach (hey, it was good enough for Maxwell!).