headers
Stavros Macrakis | 1 Jun 2007 01:49
Picon
Favicon
Gravatar

Re: jacobian not listed

On 5/31/07, Barton Willis <willisb <at> unk.edu> wrote:

Sometimes, I think that more of Maxima should be written in Maxima.

 As Wolfgang suggests, part of it is comfort with the language.  As a language, I generally prefer Lisp to Maxima, though as soon as you start writing mathematical formulae, Maxima's infix syntax makes it handier.

The speed advantages of CL over Maxima aren't all that great for symbolic work.

I think it's useful to think of Maxima as a layered system, though of course that structure isn't explicit in the code.  You have the low-level simplification, algebraic, comparison, and arithmetic routines; then the routines that use those to perform more complicated functions like equation solving, matrix operations, and integration; and finally application code which uses those to solve the user's problem.  There is also of course the Maxima programming infrastructure including both the interpreter and the translator; those are on a different dimension in some way.

I believe it is worthwhile to keep the low-level routines fast both in their algorithms and their implementation, as sometimes the scale of problem that can be solved depends on it -- not to mention simple response time for things like plotting.  It may well be useful for some symbolic calculation to assume then retract millions of assumptions over dozens of variables, for example. One component we are missing at this level is a fast, general pattern-matcher or unifier, so the middle-level routines end up doing too much of that "by hand".

In the middle routines, the quality of the algorithm tends to be more important than the details of the implementation. On the other hand, they may depend on data structures or operations (especially destructive list operations) for efficiency.  Maxima doesn't provide those, partly because they are considered too low-level and too dangerous for the non-programmer.

And finally in user code, the algorithm is pretty much all that matters.  What tends to make code slow here is intermediate expression explosion more than the micro-efficiency of the programming constructs.  For example, getting the order of variable elimination right may easily make the difference between a solvable and an unsolvable problem. Still, it would be nice if the translator were better....

That's for symbolic work.  For heavily numeric applications, things are quite different.  There are many user-written numeric problems where efficiency does matter.  And my impression is that the translator actually is tuned pretty well for those cases, as long as you include full declarations and use the right kind of arrays.

Well, those are my thoughts off the top of my head.

             -s

_______________________________________________
Maxima mailing list
Maxima <at> math.utexas.edu
http://www.math.utexas.edu/mailman/listinfo/maxima
Permalink | Reply | 0 comments |
headers
Gottfried Helms | 1 Jun 2007 18:49
Picon
Favicon

Problem with TAYLOR

Hi,

 possibly my question can be solved differently, but
 I didn't get it this way.

 With TAYLOR() I get the first six series coefficients for

  TAYLOR( (1/(1+x) , x, 0 ,6)

 This works also for m'th powers of the expression:

  TAYLOR( (1/(1+x)^m , x, 0 ,6)

 What I wanted to get is

  TAYLOR( (1/(1+x^m)^(1/m) , x, 0 ,6)

 but this gives only one coefficients and three ellipses.
 The final exponent doesn't matter, also for

  TAYLOR( (1/(1+x^m) , x, 0 ,6)

 I get only one coefficient.

How can I proceed?

Gottfried Helms
Permalink | Reply | 3 comments |
headers
Stavros Macrakis | 1 Jun 2007 21:55
Picon
Favicon
Gravatar

Re: Problem with TAYLOR

Taylor expands into a series a0 + a1*x + a2*x^2 + ...., where the exponents of x are constant integers.

It can only expand into series like a0 + a1*x^m + a2*x^(2*m) + ..., where m remains symbolic, under rather restricted circumstances.

There are two approaches you could take here.  One is to expand 1/(q+1)^(1/m) then substitute q=x^m; the other is to use the powerseries routine.

powerseries can expand some series with closed-form coefficience.  In this case, for example:

        powerseries(1/(x^m+1)^(1/m),x,0) =>
          ('sum(x^(i1*m)/beta(-1/m-i1+1,i1+1),i1,0,inf))/(1-1/m)

two-dimensional display:

    inf
    ====              i1 m
    \                x
     >     --------------------------
    /             1
    ====   beta(- - - i1 + 1, i1 + 1)
    i1 = 0        m
    ---------------------------------
                      1
                  1 - -
                      m

Using the substitution approach:

subst(x^m,q,taylor(1/(1+q)^(1/m),q,0,6)) =>
       (120*m^5+274*m^4+225*m^3+85*m^2+15*m+1)*x^(6*m)/(720*m^6)-(24*m^4+50*m^3+35*m^2+10*m+1)*x^(5*m)\
/(120*m^5)+(6*m^3+11*m^2+6*m+1)*x^(4*m)/(24*m^4)-(2*m^2+3*m+1)*x^(3*m)/(6*m^3)+(m+1)*x^(2*m)/(2*m^2)-x\
^m/m+1

Hope this helps,

            -s

On 6/1/07, Gottfried Helms <helms <at> uni-kassel.de> wrote:
Hi,

possibly my question can be solved differently, but
I didn't get it this way.

With TAYLOR() I get the first six series coefficients for

  TAYLOR( (1/(1+x) , x, 0 ,6)

This works also for m'th powers of the expression:

  TAYLOR( (1/(1+x)^m , x, 0 ,6)

What I wanted to get is

  TAYLOR( (1/(1+x^m)^(1/m) , x, 0 ,6)

but this gives only one coefficients and three ellipses.
The final exponent doesn't matter, also for

  TAYLOR( (1/(1+x^m) , x, 0 ,6)

I get only one coefficient.

_______________________________________________
Maxima mailing list
Maxima <at> math.utexas.edu
http://www.math.utexas.edu/mailman/listinfo/maxima
Permalink | Reply | 2 comments |
headers
Gottfried Helms | 2 Jun 2007 14:02
Picon
Favicon

Re: Problem with TAYLOR

Am 01.06.2007 21:55 schrieb Stavros Macrakis:
> Taylor expands into a series a0 + a1*x + a2*x^2 + ...., where the
> exponents of x are constant integers.
> 
> It can only expand into series like a0 + a1*x^m + a2*x^(2*m) + ...,
> where m remains symbolic, under rather restricted circumstances.
> 
> There are two approaches you could take here.  One is to expand
> 1/(q+1)^(1/m) then substitute q=x^m; the other is to use the powerseries
> routine.
> 
> powerseries can expand some series with closed-form coefficience.  In
> this case, for example:
> 
>         powerseries(1/(x^m+1)^(1/m),x,0) =>
>           ('sum(x^(i1*m)/beta(-1/m-i1+1,i1+1),i1,0,inf))/(1-1/m)
> 
> two-dimensional display:
> 
>     inf
>     ====              i1 m
>     \                x
>      >     --------------------------
>     /             1
>     ====   beta(- - - i1 + 1, i1 + 1)
>     i1 = 0        m
>     ---------------------------------
>                       1
>                   1 - -
>                       m
> 
Hi Stavros,

 yes, this works fine, thanks.

 But now I don't know, whether I have a mathematical or
 still a technical problem...

 If I cancel in the above sum, (taking care not to use m,
 which lead to divergences) I get a simple binomial-formula
 with the main terms [1/m : i1]*(x^m)^i1
 (where ([a:b] is binomial(a,b))

 sum(i1=0,inf,  binomial(1/m, i1)*(x^m)^i1 )

 where the first parameter in binomial(1/m,i1) is constant.

 What I was expecting was, that the first parameter would be
 varying and the second constant...

 Example:

 Assume the pascal-matrix, containing the binomial-coefficients,
P=
  1
  1  1
  1  2  1
  1  3  3  1
  ...
PInv = P^-1 =
  1
 -1  1
  1 -2  1
 -1  3 -3  1
  ...

and the powerseries-vector V(x) = column([1,x,x^2,x^3,... ])

Then the result expresses, with m=1,

   P * V(x) = V(1+x)

 and the first parameters in the binomials of a row of P are constant,
 reflecting the maxima-result.

 What I wanted was (if m=1)

 1/x* V(1/x))~ * P = 1/(1+x) * V(1/(1+x))~

 and here the coefficients of P are read column-wise, thus the
 first parameter of the binomials is varying and the second is
 constant.

 So, to define the generating-function for the first column of
 P I use (with approprate scaling)
          x
 TAYLOR( --- ,x,0,4) = x V(x)~ * column([1,-1,1,-1,...])
         1+x

 for the second column I use

           x
 TAYLOR( (---)^2,x,0,4 ) =  x^2 V(x)~ * column([1,-2,3,-4,...])
          1+x

 for the second column I use

           x
 TAYLOR( (---)^3,x,0,6 ) =  x^3 V(x)~ * column([1,-3,6,-10,...])
          1+x

 and so on, where binomials are produced, which are varying in the first
 parameter and constant in the second.

 The exponent m at x then should provide a variability of the stepwidth, and
 I thought it would be an appropriate idea to replace x by x^m and revert this
 by powering of the whole expression by 1/m.

 But how can then the form of the resulting formula change to an
 expression, which reflects a P*V(x)-multiplication instead a
 V(1/x)~ * P - multiplication???
 Or is this operation ambiguous, are both ways possible and maxima
 just gets the simpler one?

 * scratchhead *

Gottfried Helms
Permalink | Reply | 1 comment |
headers
Gottfried Helms | 2 Jun 2007 14:56
Picon
Favicon

Re: Problem with TAYLOR / corrected text

(Resend it because of some dumb textual mistakes:)

Am 01.06.2007 21:55 schrieb Stavros Macrakis:
> Taylor expands into a series a0 + a1*x + a2*x^2 + ...., where the
> exponents of x are constant integers.
>
> It can only expand into series like a0 + a1*x^m + a2*x^(2*m) + ...,
> where m remains symbolic, under rather restricted circumstances.
>
> There are two approaches you could take here.  One is to expand
> 1/(q+1)^(1/m) then substitute q=x^m; the other is to use the powerseries
> routine.
>
> powerseries can expand some series with closed-form coefficience.  In
> this case, for example:
>
>         powerseries(1/(x^m+1)^(1/m),x,0) =>
>           ('sum(x^(i1*m)/beta(-1/m-i1+1,i1+1),i1,0,inf))/(1-1/m)
>
> two-dimensional display:
>
>     inf
>     ====              i1 m
>     \                x
>      >     --------------------------
>     /             1
>     ====   beta(- - - i1 + 1, i1 + 1)
>     i1 = 0        m
>     ---------------------------------
>                       1
>                   1 - -
>                       m
>
Hi Stavros,

 yes, this works fine, thanks.

 But now I don't know, whether I have a mathematical or
 still a technical problem...

 If I cancel in the above sum, (taking care not to use m,
 which lead to divergences) I get a simple binomial-formula
 with the main terms of binomial(1/m , i1)

 as
        sum(i1=0,inf,  binomial(1/m, i1)*(x^m)^i1 )

 where the *first* parameter in binomial(1/m,i1) is constant.

 What I was expecting was, that the first parameter would be
 *varying* and the second would be constant...

 Example:

 Assume the pascal-matrix, containing the binomial-coefficients,
P=
  1
  1  1
  1  2  1
  1  3  3  1
  ...
PInv = P^-1 =
  1
 -1  1
  1 -2  1
 -1  3 -3  1
  ...

and the powerseries-vector V(x) = column([1,x,x^2,x^3,... ])

Then the result expresses, with m=1,

   P * V(x) = V(1+x)

 and the first parameters in the binomials, taken from a row of P, are constant,
 which is reflecting the maxima-result.

 What I wanted was (if m=1)

 1/x* V(1/x))~ * PInv = 1/(1+x) * V(1/(1+x))~

 and here the coefficients of PInv are read column-wise, thus the
 first parameter of the binomials is varying and the second is
 constant.

 So, to define the generating-function for the first column of
 PInv I use (with approprate scaling)

          x
 TAYLOR( --- ,x,0,4) = x V(x)~ * column([1,-1,1,-1,...])
         1+x

 for the second column I use

           x
 TAYLOR( (---)^2,x,0,4 ) =  x^2 V(x)~ * column([1,-2,3,-4,...])
          1+x

 for the second column I use

           x
 TAYLOR( (---)^3,x,0,6 ) =  x^3 V(x)~ * column([1,-3,6,-10,...])
          1+x

 and so on, where binomials are produced, which are varying in the first
 parameter and constant in the second.

 The exponent m at x then should provide a variability of the stepwidth,
 when PInv(m) is understood as an operator, which transfers 1/x V(1/x) to 1/y V(y),
 where y = x+1 or varying stepwidths like y^m = x^m+1
 and I thought it would be an appropriate idea to replace x by x^m and revert this
 by powering of the whole expression by 1/m.

 But how can then the form of the resulting formula change to an
 expression, which reflects a P*V(x)-multiplication instead a
 V(1/x)~ * PInv - multiplication???

 Or is this operation ambiguous, are both ways possible and maxima
 just gets the simpler one?

 * scratchhead *

Gottfried Helms
Permalink | Reply | 0 comments |
headers
sen1 | 3 Jun 2007 16:39
Picon

join N lists ?

Hello,
  I need to have a function which takes a list L of lists M_i  and forms a
  long list of all the elements of the elements of the M_i's.

The function

m_join(m_list):=
listify(flatten(setify(makelist(setify(m_list[i]),i,1,length(m_list)))))$

does this.  Is there a better way?  For instance, is there a built-in
function which does some sort of join on a list of lists?

TIA,
  -sen
--

-- 
  ---------------------------------------------------------------------------
  | Sheldon E. Newhouse            |    e-mail: sen1 <at> math.msu.edu           |
  | Mathematics Department         |       				   |
  | Michigan State University      | telephone: 517-355-9684                |
  | E. Lansing, MI 48824-1027 USA  |       FAX: 517-432-1562                |
  ---------------------------------------------------------------------------
Permalink | Reply | 4 comments |
headers
Stavros Macrakis | 3 Jun 2007 17:39
Picon
Favicon
Gravatar

Re: join N lists ?

Your specification isn't entirely clear.  Do you want the set of distinct elements, or the list of all elements? Also, you use the term "join" -- does that mean something other than simply the set-theoretic union?  (It means something quite different in database theory...)

If the set of distinct elements:

       xreduce('union,map('setify, l))

or somewhat less efficiently

       setify(xreduce('append, l))
    
If the list of all elements:

       xreduce('append, l)

xreduce is a very useful function (thank you, Barton!) and is very generally useful.  On the other hand, flatten has very specialized uses, and is easy to misuse:


If any of the elements of the lists are themselves sets, it includes their elements (not the set itself), e.g.

     m_join( [ [ {a,b,c} , d], [e,f] ] ) => [a,b,c,d,e,f]    NO!
            should be [  {a,b,c} , d, e, f  ]

This is because you're using the "flatten" function.  The flatten function is only really appropriate for recursive representations.  So, for example, if you represent a tree as

       mytree: node(a,node(node(b,c),d),e)$

to get all the leaves under the top node:

       flatten(mytree) => node(a,b,c,d,e)

If, on the other hand, you have a well-defined structure, e.g. a list of lists, then flatten is just an invitation to bugs such as the above.  Flatten is very rarely the right thing.

        -s

_______________________________________________
Maxima mailing list
Maxima <at> math.utexas.edu
http://www.math.utexas.edu/mailman/listinfo/maxima
Permalink | Reply | 3 comments |
headers
Gosei Furuya | 3 Jun 2007 18:08
Picon

Re: join N lists ?

hi Stavros

I would like to use flatten with outermap,
I wrote newflatten dunction.

(%i4) load("newflatten.lisp");
(%o4)                     newflatten.lisp
(%i5) nflatten([[[a,b,c],d],[e,f]],1);
(%o5)                  [[a, b, c], d, e, f]
(%i6) nflatten([[[a,b,c],d],[e,f]]); /*default level 3*/
(%o6)                   [a, b, c, d, e, f]
(%i7) nflatten([[[a,b,c],d],[e,f]],2);
(%o7)                   [a, b, c, d, e, f]
how do you think?
gosei furuya

2007/6/4, Stavros Macrakis <macrakis <at> alum.mit.edu>:
Your specification isn't entirely clear.  Do you want the set of distinct elements, or the list of all elements? Also, you use the term "join" -- does that mean something other than simply the set-theoretic union?  (It means something quite different in database theory...)

If the set of distinct elements:

       xreduce('union,map('setify, l))

or somewhat less efficiently

       setify(xreduce('append, l))
    
If the list of all elements:

       xreduce('append, l)

xreduce is a very useful function (thank you, Barton!) and is very generally useful.  On the other hand, flatten has very specialized uses, and is easy to misuse:

If any of the elements of the lists are themselves sets, it includes their elements (not the set itself), e.g.

     m_join( [ [ {a,b,c} , d], [e,f] ] ) => [a,b,c,d,e,f]    NO!
            should be [  {a,b,c} , d, e, f  ]

This is because you're using the "flatten" function.  The flatten function is only really appropriate for recursive representations.  So, for example, if you represent a tree as

       mytree: node(a,node(node(b,c),d),e)$

to get all the leaves under the top node:

       flatten(mytree) => node(a,b,c,d,e)

If, on the other hand, you have a well-defined structure, e.g. a list of lists, then flatten is just an invitation to bugs such as the above.  Flatten is very rarely the right thing.

        -s


_______________________________________________
Maxima mailing list
Maxima <at> math.utexas.edu
http://www.math.utexas.edu/mailman/listinfo/maxima


_______________________________________________
Maxima mailing list
Maxima <at> math.utexas.edu
http://www.math.utexas.edu/mailman/listinfo/maxima
Permalink | Reply | 2 comments |
headers
Stavros Macrakis | 3 Jun 2007 18:39
Picon
Favicon
Gravatar

Re: join N lists ?

On 6/3/07, Gosei Furuya <go.maxima <at> gmail.com> wrote:

I would like to use flatten with outermap,
I wrote newflatten dunction.

(%i4) load("newflatten.lisp");
(%o4)                     newflatten.lisp
(%i5) nflatten([[[a,b,c],d],[e,f]],1);
(%o5)                  [[a, b, c], d, e, f]
(%i6) nflatten([[[a,b,c],d],[e,f]]); /*default level 3*/
(%o6)                   [a, b, c, d, e, f]
(%i7) nflatten([[[a,b,c],d],[e,f]],2);
(%o7)                   [a, b, c, d, e, f]
how do you think?

I like the idea of specifying how many levels.  In particular, this means you could give a useful error message in cases where N levels aren't present, which will, I suspect, catch a lot of errors -- but I see from %o6 that you do not do this.  I would have thought that the intuition of nflatten was that you know a priori that the object consists precisely of nested lists to N levels.  Also, I wonder if nflatten works only for the "[" operator, or for other operators as well, and for mixed operators, e.g. sets of lists or lists of sets.  If so, is the result a list or a set?

              -s
_______________________________________________
Maxima mailing list
Maxima <at> math.utexas.edu
http://www.math.utexas.edu/mailman/listinfo/maxima
Permalink | Reply | 1 comment |
headers
Gosei Furuya | 3 Jun 2007 19:02
Picon

Re: join N lists ?

hi Stavros
This cords is rewritten Barton's flatten.lisp.
So lists of sets  or sets of list  are almost  correct.
And  level error  is free. but it can not treat matrix.
(%i9) nflatten([[{a,b,c},d],[e,f]],2);
(%o9)                  [{a, b, c}, d, e, f]
(%i10) nflatten([[{a,b,c},d],{e,f}],2);
(%o10)                 [{a, b, c}, d, {e, f}]
(%i11) nflatten({[{a,b,c},d],[e,f]},2);
(%o11)                {[{a, b, c}, d], [e, f]}
(%i12) nflatten({[{a,b,c},d],[e,f]},3);
(%o12)                {[{a, b, c}, d], [e, f]}
(%i13) nflatten({{{a,b,c},d},[e,f]},2);
(%o13)                  {a, b, c, d, [e, f]}

gosei furuya

(defun $nflatten (e &optional (n 5))
  (setq e (ratdisrep e))
  (cond ((or ($atom e) ($subvarp e)(or (member ($inpart e 0) (list '&^ '&=))))
         e)
        (t
         (let ((op (multiple-value-list (get-op-and-arg e))))
           (setq e (cadr op))
           (setq op (car op))
           (setq e (mapcar #'(lambda (x) (flatten-op x op n)) e))
           (setq e (reduce #'append e))
           (cond ((and (consp (car op)) (eq (caar op) 'mqapply))
                  (append op e))
                 (t
                  `(,op , <at> e)))))))

(defun flatten-op (e op nlev)
  (let ((e-op) (e-arg))
    (setq e-op (multiple-value-list (get-op-and-arg e)))
    (setq e-arg (cadr e-op))
    (setq e-op (car e-op))
    (cond ((and (>= nlev 1)(equal e-op op))
           (mapcan #'(lambda (x) (flatten-op x op (1- nlev))) e-arg))
          (t
           (list e)))))


 
 

2007/6/4, Stavros Macrakis <macrakis <at> alum.mit.edu>:
On 6/3/07, Gosei Furuya <go.maxima <at> gmail.com> wrote:
I would like to use flatten with outermap,
I wrote newflatten dunction.

(%i4) load("newflatten.lisp");
(%o4)                     newflatten.lisp
(%i5) nflatten([[[a,b,c],d],[e,f]],1);
(%o5)                  [[a, b, c], d, e, f]
(%i6) nflatten([[[a,b,c],d],[e,f]]); /*default level 3*/
(%o6)                   [a, b, c, d, e, f]
(%i7) nflatten([[[a,b,c],d],[e,f]],2);
(%o7)                   [a, b, c, d, e, f]
how do you think?

I like the idea of specifying how many levels.  In particular, this means you could give a useful error message in cases where N levels aren't present, which will, I suspect, catch a lot of errors -- but I see from %o6 that you do not do this.  I would have thought that the intuition of nflatten was that you know a priori that the object consists precisely of nested lists to N levels.  Also, I wonder if nflatten works only for the "[" operator, or for other operators as well, and for mixed operators, e.g. sets of lists or lists of sets.  If so, is the result a list or a set?

              -s
 

_______________________________________________
Maxima mailing list
Maxima <at> math.utexas.edu
http://www.math.utexas.edu/mailman/listinfo/maxima
Permalink | Reply | 0 comments |

Gmane