symbolic vector calculus
jdo...@gmail.com
I have a simple question about the capabilities of Sage that I have
not been able to resolve by looking at the documentation. I often
find myself manipulating somewhat complex functions that take vector
arguments. I then need to derive gradients, hessians, etc. I need to
do this with out knowing the dimensions of the vectors. So for
example, what I would like to do is something like the following.
First, I would define the function f(x)=.5 * x' * A * x + b, perhaps
something like:
A = matrix();
x = vector();
b = vector();
f = function( x' * A * x + b);
Then, I would like Sage to do the calculus for me, something like,
say:
f.gradient()
sage) 2*A*x
f.hessian()
sage) A
(Of course, I wouldn't bother using a computer algebra system for such
a simple function, but you get the idea.) My question is: can Sage do
things like this? I would like to avoid specifying the dimensions of
x, or giving the entries of A when I define the function.
Thanks!
Justin
William Stein
No, Sage can't do that. There is an ambitious project by Gary Furnish
to improve the symbolic manipulation capabilities of Sage, but I think
even that wouldn't at all address the problem you state above, because
you don't want to specify that dimensions.
Of course, if I wanted to do the above in Sage I would likely define
some new Python classes that implement that sort of thing (i.e., I
would just write it).
-- William
Gary Furnish
My system does not currently support this, and I have no plans to
implement this, but it could be done pretty easily if there was a
dimensionless vectorspace and matrixspace object, but I'm not
completely convinced I see the point. If one wanted to do this they
could just choose an arbitrary dimension, and as long as they don't do
anything that depends on the dimension it should still give the same
answer (in my system, not the current one).
Tim Lahey
> My system does not currently support this, and I have no plans to
> implement this, but it could be done pretty easily if there was a
> dimensionless vectorspace and matrixspace object, but I'm not
> completely convinced I see the point. If one wanted to do this they
> could just choose an arbitrary dimension, and as long as they don't do
> anything that depends on the dimension it should still give the same
> answer (in my system, not the current one).
The purpose for this is computations like this come up in Finite Element
Analysis when you are deriving the system of equations (and probably
elsewhere).
The reason for not specifying the dimensions is that you don't know the
length of the vectors in advance since it is dependent upon the number
of
elements you choose. They have a definite structure, so once the size is
specified, you can then fill in the vectors, but you don't specify the
size
in advance. The most you would do is say that A is n by n and x is n
by 1.
Maxima and Maple can't do these calculations (I've tried), but
Mathematica
can. I had to go through a lot of effort to figure out an alternative
to use
in Maple for deriving Finite Element equations, although I haven't
been able
to convert it into Sage.
Cheers,
Tim Lahey
---
Tim Lahey
PhD Candidate, Systems Design Engineering
University of Waterloo
William Stein
Well I want Sage to be able to nicely do what you want.
Any chance you could make up a slightly nontrivial example
in Sage vapor-code, which would illustrate the sort of functionality
you enjoy in Mathematica and perhaps wish Sage had?
-- William
Ondrej Certik
Do you have some FEM code in Python?
We develope sfepy (CCing sfepy-devel):
http://code.google.com/p/sfepy/
and we'd like to hook symbolic capabilities to it, so that you can do
this kind of things symbolically. We just started to use SymPy for
checking the numerical solutions symbolically and my secret plan is to
be able to just write an equation in SymPy or Sage, optionally specify
some boundary conditions and get it solved using FEM.
Do you have some pointers to the manipulations in Mathematica? Let's
implement the same in Sage using the same or similar syntax.
I don't think it is difficult.
Ondrej
William Stein
I just want to quickly mention in case the originally poster doesn't
know that (1) sympy is part of Sage, and (2) Ondrej is the lead
sympy developer.
> Do you have some pointers to the manipulations in Mathematica? Let's
> implement the same in Sage using the same or similar syntax.
> I don't think it is difficult.
>
> Ondrej
>
> >
>
--
William Stein
Associate Professor of Mathematics
University of Washington
http://wstein.org
Ondrej Certik
Anyway, I think it would be extremely cool if Sage could solve stuff
using FEM. And sfepy could perfectly be hooked in, it is pure python +
some C stuff (we use swig for historical reasons, but could switch to
cython, or just leave the generated C file in there) + scipy + numpy.
And umfpack and pysparse. But those are all either small libraries or
already in Sage.
It's still too early to talk about inclusion, we first need to create
a spkg, test it, etc. etc., also I'd wait with this until sfepy
matures a bit more, but I think it'd make Sage extremely useful for
engineering applications. Imagine downloading sage on any computer,
typing "sage", then one simple command to type in the equation and
then get it solved automatically and correctly. Wow.
Before that one would type one command to download it from optional
spkg packages.
I created an issue for that:
http://code.google.com/p/sfepy/issues/detail?id=50
Ondrej
Tim Lahey
On Jul 21, 2008, at 4:17 PM, Ondrej Certik wrote:
>
> Do you have some FEM code in Python?
>
No, at the moment, it is split between Maple and MATLAB.
I use Maple to derive the mass and stiffness element matrices
and output a MATLAB function where I do the assembly and
then the solution. I derive my equations using the
Raleigh-Ritz method whereas most FEM codes I've found
appear to want the weak form derived through Galerkin. Since
Galerkin is difficult for coupled systems, I've avoided it.
I've run into roadblocks trying to implement things in
Sage because a number of the techniques I've used in Maple
don't easily translate.
> We develope sfepy (CCing sfepy-devel):
>
> http://code.google.com/p/sfepy/
>
> and we'd like to hook symbolic capabilities to it, so that you can do
> this kind of things symbolically. We just started to use SymPy for
> checking the numerical solutions symbolically and my secret plan is to
> be able to just write an equation in SymPy or Sage, optionally specify
> some boundary conditions and get it solved using FEM.
>
> Do you have some pointers to the manipulations in Mathematica? Let's
> implement the same in Sage using the same or similar syntax.
> I don't think it is difficult.
In Mathematica, you can do:
D[Transpose[a].B.c, c]
to take the derivative with respect to the vector c. It worked in
V4 of Mathematica. As I've mentioned before, neither Maple or Maxima
can do this. The Maple developers have basically told me that it
wasn't going to be done.
I know that for a course a professor had FE codes
developed in Mathematica and both the notes and the Mathematica
codes were on his web site.
I think that they are:
http://www.colorado.edu/engineering/Aerospace/CAS/courses.d/IFEM.d/
and
http://www.colorado.edu/engineering/Aerospace/CAS/courses.d/NFEM.d/Home.html
The first is for linear FEM and the second is for non-linear FEM.
However, I can't recall if they covered exactly what I'm doing.
I'm happy to share the Maple code for the FE methods I've been
developing
since manipulations like the one above will be some time in coming
(if at all). I have two example programs in addition to the basic
package.
E-mail me off-list if you want to see it and possibly help convert it
to Sage/sympy.
I've also been looking at basic Calculus of Variations (e.g., being
able to take a variation and manipulate it). For that, I need the
ability
to examine integrands, integral bounds, and manipulate individual terms
(and parts of terms). In Maple, you can do this with op and various
integration commands.
Lastly, if the ability to work with symbolic vectors/matrices was added,
a nice thing to have would be to work with partitioned matrices and
vectors. For example:
A = [[ B, C],[D, E]]
where B, C, D, and E are general sub-matrices. This crops up all the
time
in control theory.
Harald Schilly
> Lastly, if the ability to work with symbolic vectors/matrices was added,
> a nice thing to have would be to work with partitioned matrices and
> vectors. For example:
>
> A = [[ B, C],[D, E]]
>
> where B, C, D, and E are general sub-matrices.
sage: A = matrix(QQ, 2, 2, [3,9,6,10])
sage: block_matrix([A, -A, ~A, 100*A])
[ 3 9| -3 -9]
[ 6 10| -6 -10]
[-----------+-----------]
[-5/12 3/8| 300 900]
[ 1/4 -1/8| 600 1000]
Harald
Tim Lahey
Thanks, Harald. But I was talking about in the more general sense that
B, C, D, and E are general (unknown elements) matrices at the time of
manipulation.
Many manipulations that are simple to do by hand have little to no
support in various CAS programs. Partitioned matrices happen to be one
of them. Others include: differentiation with respect to a function
(I get around that in Maple by doing a substitution with a symbol and
then reversing the substitution), derivatives with respect to a vector
(mentioned before in this thread), and Calculus of Variations (I've done
work on implementing some of it in Maple).
Ondrej Certik
Yes, those are very famous notes, I read them carefully couple years ago.
>
> I'm happy to share the Maple code for the FE methods I've been
> developing
> since manipulations like the one above will be some time in coming
> (if at all). I have two example programs in addition to the basic
> package.
> E-mail me off-list if you want to see it and possibly help convert it
> to Sage/sympy.
Unfortuntaley I am too busy to do this myself. But if you'd be
interested to do the work, I am sure you'll find a lot of enthusiastic
people willing to help you prepare the patches.
Basically, unless I am wrong, all you need are just non commutative
symbols, no?
E.g. in SymPy one can do already:
In [1]: A = Symbol("A", commutative=False)
In [2]: x = Symbol("x", commutative=False)
In [3]: x*A*x
Out[3]: x*A*x
In [4]: (x*A*x).diff(x)
Out[4]: A*x + x*A
We are missing the abstract transposition operation in sympy (we just
do it immediatelly on matrices). But that's just a matter of
implementing one more class.
Could you please post your code somewhere on the web, so that we can
all look at it and discuss it publicly? More people more opinions.
I don't know, maybe my naive approach above would not work. If you
could post here the exact operations and results that you need, it'd
be then easy to see if the above works or not.
>
> I've also been looking at basic Calculus of Variations (e.g., being
> able to take a variation and manipulate it). For that, I need the
> ability
> to examine integrands, integral bounds, and manipulate individual terms
> (and parts of terms). In Maple, you can do this with op and various
> integration commands.
>
> Lastly, if the ability to work with symbolic vectors/matrices was added,
> a nice thing to have would be to work with partitioned matrices and
> vectors. For example:
>
> A = [[ B, C],[D, E]]
>
> where B, C, D, and E are general sub-matrices. This crops up all the
> time
> in control theory.
It really depends on the particular application.
Ondrej
jdo...@gmail.com
Thank you all for the discussion. Let me give a more realistic
example of what I would like to do. Here, I use the single quote for
a transpose.
g(v) = (A'v + w)' exp(B'v + u - 1)
so, v and u are vectors, and A and B are matrices. Now the gradient
of g is (if I worked this out correctly):
dg(v)/dv = exp(B'v + u - 1) 1' . B' + exp(B'v + u - 1) 1' . B' . (B'v
+ u - 1) 1'
Here, I use " . " for the element-wise product, and " 1' " to
represent a row vector of ones. I don't believe it is possible to
represent the output with out both of these tricks.
I hope this helps to get a feeling for what would need to be done in
order to do these things automatically. It seems very challenging to
me. I actually checked out Mathematica on the recommendation of Tim
Lahey above, and while it can handle simple cases like x' A x, so far
as I can tell, it can't handle things like the above.
Justin
Ondrej Certik
It is challenging. But I would like Sage or SymPy to be able to do
that, so I created a new issue for that:
http://code.google.com/p/sympy/issues/detail?id=959
I really want the computer to do automatically anything that I am
doing by a pen and a paper.
Ondrej
William Stein
http://justindomke.wordpress.com/2008/07/22/vector-calculus/
-- William
--
Robert Dodier
> In Mathematica, you can do:
>
> D[Transpose[a].B.c, c]
>
> to take the derivative with respect to the vector c. It worked in
> V4 of Mathematica. As I've mentioned before, neither Maple or Maxima
> can do this. The Maple developers have basically told me that it
> wasn't going to be done.
If someone wants to help, maybe we can work on it together.
Part of my pitch for Maxima is that the combination of symbolic
and numerical methods is well-suited for practical problems;
FEM is pretty much a perfect example.
FWIW
Robert Dodier
Maxima developer
Sherif
I came to post a similar question and stumbled upon this discussion, which I realize is from 2008 - perhaps some of this functionality has been developed by now!
In any case, I'm brand new to Sage, as in I haven't used it at all yet though I'm excited to try it. So, if my question is answered somewhere in the documentation, you can just direct me to the right spot, but I haven't been able to find it.
Basically, I have a function of three vectors, let's call them x1, x2 and x3. The function consists of a complicated expression using dot and cross products. I may know the dimension of the vector space, so this is not the primary issue (also as far as I know, cross products are only defined for 3 or 7 dimensions). The issue is that I want the answer to be in terms of the vectors, not their individual components / coordinates in an equivalent Euclidean vector space. The function is vector-valued (the output value is of the same dimension as the input vectors), and I need to apply a composition of this function multiple times. So, the result will be overly complicated and needs to be simplified using vector identities. So, I need to obtain something like f(a4, b4, f(a3, b3, f(a2, b2, f(a1, b1, c1)))).
Can Sage do anything like this currently? If not, is there a suggested workaround with Sage, and if not, then what other software might help?
Thanks so much!
Sherif