Where there is a choice, Sage should feel like traditional mathematics
when posing mathematical problems, rather than a traditional
programming language. The following interface ideas follow from that
principle. Could this be a starting point for an SEP?
Mathematicians are expected to be able to do simple algebraic
manipulations, whereas there is no such expectation for most
programming languages. In cases where Sage's current behavior doesn't
feel like traditional mathematics, the reason is either inertia (i.e.
Sage has copied the interface of packages it delegates to, many of
which can't do symbolic manipulation, and so can't have interfaces
that feel like mathematics), or deficiencies, some only historical, in
Sage's own symbolic manipulation abilities.
Most of the proposed interfaces would put a heavier burden on Sage's
symbolic capabilities. That shouldn't stop us, though. Sage is
poised to make enormous strides on this front, and a list of desirable
interfaces will serve to focus this effort. At first pass, I will
leave out implementation details, but fast_float and pynac are
promising.
Most of these ideas are lifted from Mathematica because I know that
system well, and much of its interface is well thought out. Where we
can learn something even better from other systems, we should.
* Constraints should be expressed as equations and inequalities,
rather than simple lists of numbers. For example, in
minimize_constrained, the first example is currently
Let us maximize $x+y-50$ subject to the following constraints: $50*x
+24*y<=2400$,
$30*x+33*y<=2100$, $x>=45$, and $y>=5$.
sage: y = var('y')
sage: f = lambda p: -p[0]-p[1]+50
sage: c_1 = lambda p: p[0]-45
sage: c_2 = lambda p: p[1]-5
sage: c_3 = lambda p: -50*p[0]-24*p[1]+2400
sage: c_4 = lambda p: -30*p[0]-33*p[1]+2100
sage: a = minimize_constrained(f,[c_1,c_2,c_3,c_4],[2,3])
sage: a
(45.0, 6.25)
A better interface would allow us to translate the problem
specification almost directly into code:
sage: y = var('y')
sage: nminimize(x+y-50, constraints=[50*x+24*y<=2400, 30*x
+33*y<=2100, x>=45, y>=5], initial_values={x:2, y:3})
sage: {x:45.0, y:6.25}
Notice that with this notation, there isn't any need to distinguish
between minimize and minimize_constrained. On the other hand, there
should be a distinction between symbolic and numerical algorithms,
hence nminimize. Mathematica makes a further distinction between
local and global numerical algorithms, so the hierarchy is Minimize,
NMinimize, FindMinimum, where Minimize is symbolic, NMinimize tries to
find all the minima and sometimes gives up, and the less ambitious but
more robust FindMiminum only tries to find one minimum. Similarly,
there is Solve, NSolve, and FindRoot. The distinction between the N*
and Find* seems to be a source of frustration and confusion for
beginners. Can we do better?
* Differential equations should be expressed as equations. For
example, in ode_solver, the first example is
Consider solving the Van der pol oscillator x''(t) + ux'(t)
(x(t)^2-1)+x(t)=0 between t=0 and t= 100.
As a first order system it is
x'=y
y'=-x+uy(1-x^2)
Let us take u=10 and use initial conditions (x,y)=(1,0) and use the
runga-kutta prince-dormand
algorithm.
sage: def f_1(t,y,params):
... return[y[1],-y[0]-params[0]*y[1]*(y[0]**2-1)]
sage: def j_1(t,y,params):
... return [ [0, 1.0],[-2.0*params[0]*y[0]*y[1]-1.0,-
params[0]*(y[0]*y[0]-1.0)], [0,0] ]
sage: T=ode_solver()
sage: T.algorithm="rk8pd"
sage: T.function=f_1
sage: T.jacobian=j_1
sage:
T.ode_solve(y_0=[1,0],t_span=[0,100],params=[10],num_points=1000)
sage: T.plot_solution(filename='sage.png')
Let's have an interface that lets us write down the problem directly
sage: var('u,t')
sage: function('x',t)
sage: sol = ndsolve([(diff(x,
2)+u*diff(x)*(x^2-1)+x==0).substitute(u=10),x(0) == 0, diff(x)(0) ==
1], x, (0,100), absolute_error=10^-6, algorithm="rk8pd")
sage: plot(sol[x])
There's a lot going on here, so I'll unpack things roughly in the
order they appear in the example. Above is the sparsest way I can
think to write things down, and arguably, more explicit hints should
be required. I think it's enough to uniquely define the problem,
though.
* function('x',t) doesn't currently work the way I would expect--you
need to write x = function('x',t)--but if var('u') can inject u into
the name space, then function('x',t) should inject x into the name
space in the same way for consistency and convenience.
* Please don't make me state the problem in terms of first order
equations only. Turning this second order equation into a coupled set
of first order equations is a simple (but boring) algebraic
manipulation, and Sage should be able to do it. In cases where Sage
can't pull off the necessary manipulation, fail with a descriptive
error message.
* With pynac, we'll have the substitute command, so we may as well
make use of it rather than having an explicit params keyword. It
would be nice to have a way to substitute for symbolic variables in a
block. For example:
sage: @substituting(u=10)
... sol = ndsolve([diff(x,2)+u*diff(x)*(x^2-1)+x==0, x(0) ==
0, diff(x)(0) == 1], x, (0,100), absolute_error=10.0^-6,
algorithm="rk8pd")
sage: u
u
This is like the With command in Mathematica. I may be
misunderstanding how decorators work. If so, alternatives?
* Like constraints, boundary conditions should be stated as
equations rather than lists of numbers. I'm agnostic as to whether
these belong in a list with the primary equations, or whether they
should be passed as a separate argument.
* Since we've defined x to be a function of t, we don't need to be
explicit about writing x(t) or t anywhere.
* The range specifications here should be consistent with those for
plot. The same holds with other functions like integral, or a
hypothetical nintegral. Compare to
sage: x(t) = 1/2*(1 + cos(t))
sage: plot(x, (-pi,pi))
Maybe plot is a little too permissive in this regard, but whatever
we do, it should be consistent.
* I shouldn't have to specify the number of points to use in the
solution unless I really want to. GSL and other solvers know how to
choose steps in order to try to satisfy a given error requirement.
Choose a smart default here, and allow me to ignore all that on a
first pass.
* Since the equation to be solved is symbolic, Sage should know how
to calculate the Jacobian and Hessian and pass them along to the
solver.
* Let's make the value returned by ndsolve really useful. I think
this is an area where Mathematica falls down. It returns a nested
list of replacement rules, each pointing to an interpolating
function. I know why it does this, but it's not immediately obvious
how to plot that, or even how to evaluate it at a particular value.
What you have to do is something like Plot[x[t] /. First[sol],{t,
0,100}]. If you want it to be fast, you have to do
Plot[Evaluate[x[t] /. First[sol]],{t,0,100}]. Ick. The best I can
think of is to have the return value be an object that responds to the
dictionary interface, but also includes all kinds of information about
the solution procedure accessible as properties or methods. sol[x] is
an interpolation, so I can do sol[x](pi), for example and it will spit
out a value. This interpolation object should know how to plot
itself, including by default plotting itself over its entire range,
hence plot(sol[x]). I'm pleasantly surprised to see that
SymbolicFunctionEvaluation objects can already be used as dictionary
keys, so that, for example, I can already do
sage: x = function('x',t)
sage: dict = {x:'Sweet, this actually works!'}
sage: dict[x]
'Sweet, this actually works!'
What other useful things belong in this return object?
sage: sol.function_evaluations
820
sage: sol.input_equations
[diff(x,2)+10*diff(x)*(x^2-1)+x==0, x(0)==0, diff(x)(0)==1]
sage: sol.manipulated_equations
[diff(xm[0])==xm[1], diff(xm[1])==-xm[0]+10*xm[1]*(1-xm[0]^2),
xm[0](0)==0, xm[1](0)==1]
In the above, xm is based on x, mangled in whatever way necessary
so as not to conflict with any other variables in our namespace.
Since I propose doing automatic manipulations of the equations to put
them in the required form, it would be very handy to be able to do a
reality check on these automatic manipulations.
A more explicit alternative syntax would be something like
sage: var('u,t')
sage: x = function('x')
sage: sol = ndsolve([(diff(x(t),t,
2)+u*diff(x(t),t)*(x(t)^2-1)+x(t)==0).substitute(u=10),x(0) == 0,
diff(x(t),t)(0) == 1], x(t), (t,0,100), absolute_error=10^-6,
algorithm="rk8pd")
This has the advantage of letting one write a retarded equation:
sage: sol = ndsolve([diff(x(t),t)==-x(t-2), x(0)==1], x(t), (t,0,5))
although I don't know any package that can actually solve such
systems.
Whatever we do, let's try to leave room for a future implementation of
partial differential equations.
A small drawback of the proposed interface is that it encourages
writing a ton of stuff on one line. Maybe it would be nice to be able
to do something like
sage: nds = ndsolver()
sage: nds.equations = (diff(x,
2)+u*diff(x)*(x^2-1)+x==0).substitute(u=10)
sage: nds.boundary_conditions = [x(0) == 0, diff(x)(0) == 1]
sage: #etc.
sage: sol = nds.solve()
although if we do this here, we should do it everywhere.
Additionally, I'm mildly allergic to having the explicit nds.solve()
call. If you like the way this looks, a simple work around would just
be to build a dictionary of the necessary arguments and then pass it
in when you're done.
* When Sage doesn't know the answer, it should say so in the most
useful possible way, rather than pretend there is no answer. A recent
example from the mailing list:
sage: var('x,y')
sage: solve([x==0, 1-exp(y)==0],x,y)
[]
I know this particular example is more a bug/misfeature in Maxima than
a design decision in Sage, but I think the general point is valid.
Are there better examples of places Sage can improve in this regard?
* Sage should have a pleasant interface for stating assumptions. I
don't have a lot of good ideas here, though.
I hope all this serves as a good conversation starter.
Regards,
JM
--~--~---------~--~----~------------~-------~--~----~
To post to this group, send email to sage-devel@googlegroups.com
To unsubscribe from this group, send email to [EMAIL PROTECTED]
For more options, visit this group at http://groups.google.com/group/sage-devel
URLs: http://www.sagemath.org
-~----------~----~----~----~------~----~------~--~---
