mentat-collective / emmy

The Emmy Computer Algebra System.
https://emmy.mentat.org
GNU General Public License v3.0
405 stars 24 forks source link

Question arising from vector calculus #173

Closed alexgian closed 5 months ago

alexgian commented 6 months ago

I know that Emmy already has vector calculus functions, including div and curl, but those definitions are in a DG context, which I do not yet quite follow. It might be nice to have them wrapped in a more basic form (with sensible defaults for a conventional ℝ³ user) who may want to try out some vector calc.

However, it also occurred to me as I was thinking about it, that the definition of curl, ∇✕ , could be made to match the common textbook one fairly easily:

(defn curl [F]
   (determinant 
      (matrix-by-row
         [(up 1 0 0) (up 0 1 0) (up 0 0 1)]
         [(∂ 0)  (∂ 1)  (∂ 2)]
         [(ref F 0)  (ref F 1)  (ref F 2)])))

This could be implemented by Emmy, if a couple of small changes were made.
First, redefine multiplication of an operator by a function as the application of that operator to the function. Currently it seems to work like this: (((* D sin) 'x) 'y) is interpreted as (* 'x ((D sin) 'y)) i.e (* 'x (cos 'y)) semantics which aren't clear to me, but certainly do follow the behaviour of scmutils Personally, I'd have thought that (* D sin) would be equivalent to (D sin), i.e. a simple application. Of course, there may be some important reasons that it is the way it is, but I can't see them. If anyone can explain....

Secondly, we'd also need to define multiplication of a tuple of reals by a function as the scaling of that function into the tuple, i.e. (* (up 0 1 2) sin) => (up 0 sin (* 2 sin)) At the moment this is not a case that's catered for, but I think it might make sense. With these two small changes the common determinant definition of curl is quite feasible. Just musing...

sritchie commented 5 months ago

I'll address these in two comments!

In the current implementation, operators can combine with non-operators ("co-operators") using +, - and *; I don't think I have the rules for this written down anywhere but I was careful to stick to the scmutils conventions and I don't think we can change them without getting into trouble mathematically.

Quick overview:

+ and -

When you add (or subtract) two operators, just like when you add two functions, you get a new operator that applies both operators and then adds (or subtracts) the two results:

emmy.env> (((+ D D) sin) 'x)
(* 2 (cos x))

When you add or subtract an operator with a non-operator, the non-operator is converted to an operator that composes the non-operator with the new operator's argument. (These two are covered by combine-op-f and combine-f-op in emmy.operator, and the docstring describes what happens: https://github.com/mentat-collective/emmy/blob/main/src/emmy/operator.cljc#L229-L277)

emmy.env> (((+ D cos) sin) 'x)
(+
   ;; this first part is ((D sin) 'x)
   (cos x)
   ;; the `cos` composed with `sin`, then the result applied to `x`
   (cos (sin x)))

*

Multiplication of two operators is the same as composition. This is important, so that ((square D) f) == (D (D f)).

Multiplication of an operator by a non-operator works differently depending on if the operator is on the left or the right.

Op on the left

If the operator is on the left, Emmy handles (* <operator> <co-operator>) by returning a new operator that does the following internally:

(fn [x]
  (<operator> (* <co-operator> x)))

So your (* D sin) example above becomes the operator (fn [x] (D (* sin x))), and the full example unfolds like this:

(((* D sin) 'x) 'y)
(((fn [x] (D (* sin x))) 'x) 'y)
((D (* sin 'x)) 'y)

;; remember now that (* sin 'x) is the same as `(fn [arg] (* (sin arg) ((constantly 'x) arg)))`

;; final result:
(* x (cos y))

Operator on the right

If the operator is on the right, THEN things work more like what you predicted, where (* <co-operator> <operator>) becomes

(fn [x] (* <co-operator> (<operator> x)))

So this would fail with flipped args, because (D 'x) is meaningless:

emmy.env> (((* sin D) 'x) 'y) 
Execution error (IllegalArgumentException) at emmy.calculus.derivative/fn (derivative.cljc:459).
No method in multimethod 'partial-derivative' for dispatch value: [clojure.lang.Symbol clojure.lang.PersistentVector]
sritchie commented 5 months ago

I tried for a while to get your version of the curl to work, but especially after reading around (e.g. https://math.stackexchange.com/questions/3999040/is-curl-actually-a-cross-product) I'm not convinced that the cross-product definition of curl can work... I don't think multiplication can make sense as application.

I do have this in emmy.calculus.vector-calculus:

(def Curl
  "Operator that takes a function `f` and returns a function that
  calculates the [Curl](https://en.wikipedia.org/wiki/Curl_(mathematics)) of `f`
  at its input point.

  `f` must be a function from $\\mathbb{R}^3 \\to \\mathbb{R}^3$."
  (-> (fn [f-triple]
        (let [[Dx Dy Dz] (map d/partial [0 1 2])
              fx (ref f-triple 0)
              fy (ref f-triple 1)
              fz (ref f-triple 2)]
          (up (- (Dy fz) (Dz fy))
              (- (Dz fx) (Dx fz))
              (- (Dx fy) (Dy fx)))))
      (o/make-operator 'Curl)))

which is equivalent to what you've got, just without the explicit cross-product written out.

Who knows, maybe Sussman has thoughts! Send him an email at gjs@mit.edu, the question's just as applicable to scmutils and I'm curious to see what he thinks.

sritchie commented 5 months ago

I'll close now but feel free to re-open, or continue the discussion on the closed issue...

alexgian commented 5 months ago

Thanks Sam

The existence of your Curl function covers my immediate need, so I'm happy to leave the issue closed for now.

As for the rest of it, i.e. whether the cross-product notation is mathematically valid or just a useful hand-waving mnemonic, or simply only applicable to coords in ℝ³, I'm going to have to read up a bit more before I grasp completely.
Having read that math.stackexchange q, I suspect the latter.

Thanks for the pointers.

alexgian commented 5 months ago

As for the (* <operator> <co-operator>) explanation, OK, I can accept that ((* D f) g) => (D (* f g)) rather than ((D f) g), but I'll have to meditate a bit more to be happy with it.