D[] notation input for ODEs

[john.hoebing]

After a day spent trying to implement a workaround for a special case
of this issue:

http://trac.sagemath.org/sage_trac/ticket/6480

I thought it would be a good time to review the status of some of the
earlier issues with how we represent ordinary (partial) derivatives.
Most relevant is this old thread:

http://groups.google.com/group/sage-devel/browse_thread/thread/6074f9d1b04b39f8/

which highlights the lack of a practical inverse for the D[]
notation. For instance, in sage:

sage: a=diff(f,x,x)+diff(f,x)/x
sage: str(a)
'D[0](f)(x, y)/x + D[0, 0](f)(x, y)'

but as far as I could find by reading the trac, groups, asksage and
google, there is no way to _enter_ a differential expression from such
a str() output. That capability would be important in, for example, a
javascript JSON structure for communicating an O(P)DE, or for input
from function generators in other languages, or anyplace where its
inconvenient to generate a correct diff() python/sage expression. If
I missed such a facility, my apologies, but if not, has there been any
progress on any of these issues? (Also, I wasn't very successful in
finding good JSON equation formats via google or stackoverflow, and
I'd really like to know if there are any symbolic O(P)DE formats for
javascript.)

John Hoebing

[Nils Bruin]

On Mar 31, 5:14 pm, "john.hoebing" <jlhs...@gmail.com> wrote:
> sage: a=diff(f,x,x)+diff(f,x)/x
> sage: str(a)
> 'D[0](f)(x, y)/x + D[0, 0](f)(x, y)'

If I understand correctly, you would like to be able to put the above
string into sage and get the expression back? That is of course a very
reasonable goal. The class FDerivativeOperator can be used to create
the appropriate expressions, but suffers from the fact it uses a
different syntax. We need to build a wrapper that does allow the
indexing syntax above. Something along these lines would do the trick:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
from sage.symbolic.operators import FDerivativeOperator
class Doperator:
def __init__(self,vars=None):
self.vars= [] if vars is None else vars

def __call__(self,f):
return FDerivativeOperator(f,self.vars)

def __getitem__(self,i):
if isinstance(i,tuple):
newvars=self.vars+list(i)
else:
newvars=self.vars+[i]
return Doperator(newvars)

D=Doperator()
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

With that code in place we can indeed do:

sage: var('x,y')
sage: D[0](f)(x, y)/x + D[0, 1](f)(x, y)
D[0](f)(x, y)/x + D[0, 1](f)(x, y)

Perhaps this should go into the library somewhere. I'm not sure we can
afford to predefine D like this, given it's such a commonly used
symbol. Maple does it, though.

[Harald Schilly]

On Sunday, April 1, 2012 3:15:14 AM UTC+2, Nils Bruin wrote:

sage: D[0](f)(x, y)/x + D[0, 1](f)(x, y)

For me, this just looks quite awkward. What about 

D(f, 0) or D(f, 0, 1) ?

which can also be like f.D(0) and f.D(0,1) ?

h 

[Nils Bruin]

On Apr 1, 4:56 am, Harald Schilly <harald.schi...@gmail.com> wrote:
> > sage: D[0](f)(x, y)/x + D[0, 1](f)(x, y)
>
> For me, this just looks quite awkward.

OK, some reasons *for* the notation then.

1) It's how sage prints the expression. It's very desirable to ensure
that output is also valid input.

2) If you read the square brackets as "subscript" (not an unusual
convention) then it is a direct reflection of how differential
operators are often written: D_x, D_{xx} etc. Here we have to settle
for positional indicators, because the parameter positions can't be
uniquely identified by a variable name (compare D_x(f)(x,y) and D_x(f)
(y,x) with D[0](f)(x,y) and D[0](f)(y,x) )

3) D[0] really is something mathematical: It's a differential
operator. You can't express these directly with the other notations. I
guess you could do

D0 = lambda f: diff(f(x),x).function(x)

but then sage would know very little of the nature of D0 (not to
mention that it changes what arguments of the function are)

4) [I'd say this one would need more work] If we allow arithmetic with
these you could write

Lambda = D[0,0]+D[1,1]
Dt = D[2]

(Lambda - Dt)(f) == 0

to express the heat equation. I think this would be better implemented
in a stricter, more well-defined way as
"the differential operator ring on differentiable functions over K in
variables (x,y,t)" or something like that, but that might be too
algebraic an approach for the people who would be using it most.

> What about D(f, 0) or D(f, 0, 1) ?

That's essentially the syntax we already have for FDerivativeOperator.

> which can also be like f.D(0) and f.D(0,1) ?

I think that has merit, but it doesn't solve the problem that sage
should probably have some notational convention for the mathematical
entities called differential operators.

[Michael Orlitzky]

On 04/01/2012 02:06 PM, Nils Bruin wrote:
>
> That's essentially the syntax we already have for FDerivativeOperator.
>
>> which can also be like f.D(0) and f.D(0,1) ?
>
> I think that has merit, but it doesn't solve the problem that sage
> should probably have some notational convention for the mathematical
> entities called differential operators.
>

If you think you want derivatives right now, you're probably better off
using aptly-named symbolic variables, and performing the differentiation
later on once you're done manipulating the "derivatives."

sage: f = function('f', x)
sage: f_prime = f.diff(x)
sage: f_prime.simplify()
D[0](f)(x)
sage: f_prime(0).simplify()
...
NotImplementedError: arguments must be distinct variables

Most symbolic stuff blows up the same way:

sage: solve([f_prime(0) == 0], f_prime(0))
...
NotImplementedError: arguments must be distinct variables

Substitution doesn't even work pre-evaluation:

sage: g = function('g', x)
sage: f_prime.substitute_function(f,g)
D[0](f)(x)

A symbolic var('D0fx') on the other hand will work just fine.

[john.hoebing]

On Apr 1, 11:06 am, Nils Bruin <nbr...@sfu.ca> wrote:
> On Apr 1, 4:56 am, Harald Schilly <harald.schi...@gmail.com> wrote:
>
> > > sage: D[0](f)(x, y)/x + D[0, 1](f)(x, y)
>
> > For me, this just looks quite awkward.
>
> OK, some reasons *for* the notation then.
>
> 1) It's how sage prints the expression. It's very desirable to ensure
> that output is also valid input.
>

Yes.

I agree completely about the usefulness of arithmetic on differential
operators. Another necessary feature is some way to extract
coefficients of the various D operators. For example the following
breaks with 'diff' or 'D' queries:

sage: a=f.diff(x,2)*x*x
sage: a
x^2*D[0, 0](f)(x)
sage: a.coeffs(x)
[[D[0, 0](f)(x), 2]]
sage: a.coeffs(diff(f,x))
/home/jlh/wrk/sage-5.0.beta11/local/lib/python2.7/site-packages/
IPython/iplib.py:2260: DeprecationWarning: Substitution using function-
call syntax and unnamed arguments is deprecated and will be removed
from a future release of Sage; you can use named arguments instead,
like EXPR(x=..., y=...)
exec code_obj in self.user_global_ns, self.user_ns
---------------------------------------------------------------------------
ValueError Traceback (most recent call
last)
...
ValueError: The name "D[0](f)(x)" is not a valid Python identifier.
sage:

[Nils Bruin]

On Apr 1, 2:34 pm, Michael Orlitzky <mich...@orlitzky.com> wrote:
> If you think you want derivatives right now, you're probably better off
> using aptly-named symbolic variables, and performing the differentiation
> later on once you're done manipulating the "derivatives."
>
>    sage: f = function('f', x)
>    sage: f_prime = f.diff(x)
>    sage: f_prime.simplify()
>    D[0](f)(x)
>    sage: f_prime(0).simplify()
>    ...
>    NotImplementedError: arguments must be distinct variables

Most of these can be fixed by using a set of temporary intermediate
variables. See

http://trac.sagemath.org/sage_trac/ticket/7377#comment:54

for some ideas. The restrictions that currently exist on evaluating
FDerivativeOperator expressions are largely unnecessary. Probably the
person who implemented it did the minimum amount of work to transition
from "named variable derivatives" to "position-based derivatives".

Removing these restrictions is probably a nice student project (you'll
get to see the internals of quite some parts of sage, but I don't
think it will ever be overly complicated)

[Nils Bruin]

On Apr 2, 12:25 am, Nils Bruin <nbr...@sfu.ca> wrote:

> Removing these restrictions is probably a nice student project (you'll
> get to see the internals of quite some parts of sage, but I don't
> think it will ever be overly complicated)

Turns out there is a nasty complication. Any experts on

symbolic_expression_from_maxima_string

can probably weigh in. In short, the problem is that

at( f(x,y,z), [x=u,y=v,z=w])

is valid in maxima but cannot be handled because of the list embedded
in an expression.

See http://trac.sagemath.org/sage_trac/ticket/12796

The binary-level max_to_sr and sr_to_max are easily patched to handle
these.

[Nils Bruin]

On Apr 2, 8:56 pm, Nils Bruin <nbr...@sfu.ca> wrote:
> In short, the problem is that
>
> at( f(x,y,z), [x=u,y=v,z=w])
>
> is valid in maxima but cannot be handled because of the list embedded
> in an expression.

Problem solved by extending the grammar accepted by
sage.misc.parser.Parser. Now an "arg" can also be a list.

Patch up for review on:

http://trac.sagemath.org/sage_trac/ticket/12796

[Nils Bruin]

On Apr 1, 2:34 pm, Michael Orlitzky <mich...@orlitzky.com> wrote:

> Substitution doesn't even work pre-evaluation:
>

> sage: f = function('f', x)
> sage: f_prime = f.diff(x)

>    sage: g = function('g', x)
>    sage: f_prime.substitute_function(f,g)
>    D[0](f)(x)

That's because you're trying to substitute expressions, not functions:

sage: f,g
(f(x), g(x))

For that, you would have to use a different substitution function and
it wouldn't have any effect anyway because the expression f(x) does
not occur in f_prime.

If you fix this, the substitution does work:

sage: f = f.operator()
sage: g = g.operator()
sage: f,g
(f, g)
sage: f_prime.substitute_function(f,g)
D[0](g)(x)

There is a different issue, which is fixed by http://trac.sagemath.org/sage_trac/ticket/12801

[Michael Orlitzky]

On 04/04/12 15:39, Nils Bruin wrote:
> On Apr 1, 2:34 pm, Michael Orlitzky <mich...@orlitzky.com> wrote:
>
>> Substitution doesn't even work pre-evaluation:
>>
>> sage: f = function('f', x)
>> sage: f_prime = f.diff(x)
>> sage: g = function('g', x)
>> sage: f_prime.substitute_function(f,g)
>> D[0](f)(x)
>
> That's because you're trying to substitute expressions, not functions:
>
> sage: f,g
> (f(x), g(x))
>

I sort of buy that, but if,

sage: f = function('f', x)

doesn't make `f` a function, that's a user-interface WTF =)

It wouldn't matter if I could substitute them as expressions, but that
doesn't work either:

sage: f = function('f', x)

sage: g = function('g', x)

sage: f.diff(x).substitute_expression(f==g)
D[0](f)(x)
sage: f.diff(x).subs(f=g)
D[0](f)(x)
sage: f.diff(x)(f=g)
D[0](f)(x)

I reported a lot of these at,

http://trac.sagemath.org/sage_trac/ticket/11842

and later duped it to,

http://trac.sagemath.org/sage_trac/ticket/6480

which seems to be the first report of the problem.

> For that, you would have to use a different substitution function and
> it wouldn't have any effect anyway because the expression f(x) does
> not occur in f_prime.
>
> If you fix this, the substitution does work:
>
> sage: f = f.operator()
> sage: g = g.operator()
> sage: f,g
> (f, g)
> sage: f_prime.substitute_function(f,g)
> D[0](g)(x)
>
> There is a different issue, which is fixed by http://trac.sagemath.org/sage_trac/ticket/12801
>

After #12801, I'm having trouble reconciling these two examples:

sage: f = function('f', x)

sage: g = function('g')
sage: f.diff(x).substitute_function(f,g)
D[0](f)(x)

versus,

sage: f = function('f')
sage: g = function('g')
sage: f(x).diff(x).substitute_function(f,g)
D[0](g)(x)

It would be slightly better I think if the first example worked, but I
can live with using the second (although I never would have discovered
it on my own).

What I would /really/ like to be able to do is,

midpoint = (1/2)*( f(a) + f(b) )

and then approximate multiple functions by swapping out the symbolic `f`
for a real function like sine.

[Nils Bruin]

On Apr 6, 12:03 pm, Michael Orlitzky <mich...@orlitzky.com> wrote:
> I sort of buy that, but if,
>
>   sage: f = function('f', x)
>
> doesn't make `f` a function, that's a user-interface WTF =)

I would tend to agree, but it's a fact of life that function('f',x)
returns the symbolic function f evaluated at x. The docstring of
"function" doesn't help either. The confusion between the nature of f
and f(x) is a usual one in calculus. The vagueness actually lies more
in the nature of 'x' than the nature of 'f': Is 'x' a polynomial
variable? an arbitrary real number? an arbitrary complex number? The
identity element in the ring of differentiable functions on R?

> After #12801, I'm having trouble reconciling these two examples:
>
>   sage: f = function('f', x)
>   sage: g = function('g')
>   sage: f.diff(x).substitute_function(f,g)
>   D[0](f)(x)
>
> versus,
>
>   sage: f = function('f')
>   sage: g = function('g')
>   sage: f(x).diff(x).substitute_function(f,g)
>   D[0](g)(x)
>
> It would be slightly better I think if the first example worked, but I
> can live with using the second (although I never would have discovered
> it on my own).

Again, the result of "function('f',x)" is the expression "f(x)", which
does not occur in "D[0](f)(x)". Really, substitute_function should
throw an error if it is asked to substitute something that cannot be
viewed as a "function". Note that function('f',x)(3) throws a
deprecation warning. Life would be much clearer if that were to
finally become a proper error.

> What I would /really/ like to be able to do is,
>
>   midpoint = (1/2)*( f(a) + f(b) )
>
> and then approximate multiple functions by swapping out the symbolic `f`
> for a real function like sine.

which you *can* provided you avoid the pitfall of "function('f',x)":

sage: f=function('f')
sage: var('a,b')
(a, b)
sage: midpoint = (1/2)*( f(a) + f(b) )
sage: midpoint
1/2*f(a) + 1/2*f(b)
sage: midpoint.substitute_function(f,sin)
1/2*sin(a) + 1/2*sin(b)

You can even make that into a function of a,b
sage: mp
(a, b) |--> 1/2*f(a) + 1/2*f(b)

although substitute_function silently changes the nature of the object
[I think that's a bug]:
sage: mp.substitute_function(f,sin)
1/2*sin(a) + 1/2*sin(b)

[now we just have an expression in a,b again]

That said, I think your confusion shows that it's probably better if
function('f',x) were to be deprecated in favour of function('f')(x).
The relevant bit of information to glean from the notation is that f
takes one argument. You can already indicate the valid number of
arguments:

sage: f= function('f',nargs=3)
sage: f(1,2)
TypeError: Symbolic function f takes exactly 3 arguments (2 given)

The attribute f.number_of_arguments() doesn't seem to get much use
though:

sage: sage.symbolic.operators.FDerivativeOperator(f,[3])
D[3](f)

is happily produced but can never be successfully evaluated.

[Michael Orlitzky]

On 04/06/12 16:12, Nils Bruin wrote:
>
> which you *can* provided you avoid the pitfall of "function('f',x)":
>
> sage: f=function('f')
> sage: var('a,b')
> (a, b)
> sage: midpoint = (1/2)*( f(a) + f(b) )
> sage: midpoint
> 1/2*f(a) + 1/2*f(b)
> sage: midpoint.substitute_function(f,sin)
> 1/2*sin(a) + 1/2*sin(b)
>
> You can even make that into a function of a,b
> sage: mp
> (a, b) |--> 1/2*f(a) + 1/2*f(b)

That was a too-simple example. You can't create e.g. a cubic spline
because of the evaluated derivatives. In general the form over [-1,1]
would look like,

s(f;x) = a(x)*f(-1) + b(x)*f'(-1) + c(x)*f(1) + d(x)*f'(1)

Swapping out `f` after evaluating it at the endpoints is what causes the
biggest problems.

If you make `s` a function that takes `f` as an argument, it works, but
I need to be able to swap out `f` later for two reasons:

1) With higher order approximations, optimal splines can take a long
time to compute.

2) To calculate error bounds, I need to know the coeffients of f,
f.diff(),... What's the coefficient of sin(pi)?

>
> That said, I think your confusion shows that it's probably better if
> function('f',x) were to be deprecated in favour of function('f')(x).
> The relevant bit of information to glean from the notation is that f
> takes one argument. You can already indicate the valid number of
> arguments:
>

This does work better, I'll use it from now on, thanks.

[Nils Bruin]

On Apr 6, 2:15 pm, Michael Orlitzky <mich...@orlitzky.com> wrote:
> That was a too-simple example. You can't create e.g. a cubic spline
> because of the evaluated derivatives. In general the form over [-1,1]
> would look like,
>
>   s(f;x) = a(x)*f(-1) + b(x)*f'(-1) + c(x)*f(1) + d(x)*f'(1)
>

Doesn't this do what you want? (probably needs #12796)

sage: a=function('a',nargs=1)
sage: b=function('b',nargs=1)
sage: c=function('c',nargs=1)
sage: d=function('d',nargs=1)
sage: f=function('f',nargs=1)
sage: from sage.symbolic.operators import FDerivativeOperator
sage: fprime=FDerivativeOperator(f,[0])
sage: s=a(x)*f(-1) + b(x)*fprime(-1) + c(x)*f(1) + d(x)*fprime(1)
sage: s
a(x)*f(-1) + b(x)*D[0](f)(-1) + c(x)*f(1) + d(x)*D[0](f)(1)
sage: s.substitute_function(f,sin)
-sin(1)*a(x) + sin(1)*c(x) + cos(1)*b(x) + cos(1)*d(x)
sage: s.coeff(f(1))
c(x)
sage: s.coeff(sin(pi))
0

I suspect that you're reaching the limits of usefulness of the
symbolic ring by now, though. It's probably easier to keep things
sorted yourself rather than mashing it all together in a big symbolic
expression and hope you can fish out the useful bits later on again.
Probably throwing some linear algebra at your problem will keep the
structure much cleaner.

[Michael Orlitzky]

On 04/06/2012 10:51 PM, Nils Bruin wrote:
> On Apr 6, 2:15 pm, Michael Orlitzky <mich...@orlitzky.com> wrote:
>> That was a too-simple example. You can't create e.g. a cubic spline
>> because of the evaluated derivatives. In general the form over [-1,1]
>> would look like,
>>
>> s(f;x) = a(x)*f(-1) + b(x)*f'(-1) + c(x)*f(1) + d(x)*f'(1)
>>
>
> Doesn't this do what you want? (probably needs #12796)
>

Yes, I didn't realize that #12796 fixed this the first time I looked at
it. I've reviewed it, thanks!