On Mon, Jul 21, 2008 at 10:20 PM, William Stein <wst...@gmail.com> wrote:

On Mon, Jul 21, 2008 at 10:17 PM, Ondrej Certik <ond...@certik.cz> wrote:

...

On Mon, Jul 21, 2008 at 10:12 PM, Tim Lahey <tim....@gmail.com> wrote:

...

On Jul 21, 2008, at 3:41 PM, Gary Furnish wrote:

...

On Mon, Jul 21, 2008 at 10:00 AM, William Stein <wst...@gmail.com> wrote:

...

On Mon, Jul 21, 2008 at 6:28 PM, <jdo...@gmail.com> wrote:

...

Hello all,

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.

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.

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.

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.

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