Carl Love

Carl Love

21198 Reputation

24 Badges

8 years, 246 days
Natick, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity

These are answers submitted by Carl Love

Maple 2021 has introduced the reduce option to seq. Replacing add with seq[reduce= `+`] is very efficient, and it avoids the kernel crash caused by add. In the code below, I have made your a procedure, like Tom did. But using the reduce option means we don't need a list or to use elementwise operation (B~). This is thus more efficient (especially memorywise). (And add would be efficient also, but there's a bug at the moment.)

I've simplified and generalized your expression such that AB, and any higher order such expression can be handled by the new procedure below. As you can see, the time to evaluate over all tuples and add them is barely measurable.

Iterator:-CartesianProduct([1,1]): #Force compilation

P:= 2..4: #suffixes of m and c variables
N:= [$0..4]: #evaluation values of m variables

B:= subs(
    {_C= [$P], _V= [c||P]}, 
    local r:= add(M*~_C), s:= 1+r, t:= s-add(M);
    end proc
    seq[reduce= `+`](B(v), v= Iterator:-CartesianProduct(N$(rhs(P)-1)))
memory used=0.69MiB, alloc change=0 bytes, 
cpu time=0ns, real time=3.00ms, gc time=0ns


Tables are one of the fundamental data structures built directly into the Maple kernel. They are extremely efficient and extremely flexible. They are quite suitable as the foundation for building numerous other data structures. Indeed, they could easily be used to implement MultiSets (and perhaps they are so used).

Is this suitable?:

    proc(a::Vector, N::posint)::hfloat; 
    local i::posint; 
        Sum(a[i], i= 1..N) 
    end proc,
    resultname= "sum"
function sumreturn = sum(a, N)
  r = 0;
  for i = 1:N
    r = r + a(i);
  sumreturn = r;

All the ::... declarations could be omitted. An arrow procedure could be used:

    (a, N)-> Sum(a[i], i= 1..N), 
    resultname= "sum"
Warning, procedure/module options ignored
function sumreturn = sum(a, N)
  r = 0;
  for i = 1:N
    r = r + a(i);
  sumreturn = r;

The warning can be ignored in this case.

What you're trying to do is trivial with Iterator:-CartesianProduct:

P:= 4: M:= [m||(2..P)]: N:= 4:
    for c in Iterator:-CartesianProduct([$0..N]$(P-1)) do 
        eval(B, M=~[seq](c)) 


sum2N:= proc(N::nonnegint) 
local summation:= 0, i; 
    for i from 0 to N do summation:= summation+(N+i)^2 end do
end proc;

The technique that I'm about to provide will work for any numeric dsolve solution, regardless of the number of equations, the order of the derivatives, or whether it's IVP, BVP, or DAE. First, the ODEs need a slight rearrangement. In your case, that's

deq1:= diff(y(t), t, t)= -y(t)*abs(y(t));

All that I've done is solve the ODE for its highest-order derivative. Once that is done, any derivative of y of any integer order j (including j=0 or j=1), can be expressed as 

eval[recurse](diff(y(t), [t$j]), {deq1})

and these expressions will be recognized by odeplot. The only modification to the code that is needed is to solve the ODEs for their highest-order derivatives. So, for example,  you can make an array plot of y and its first 4 derivatives by

deq1:= diff(y(t), t, t)= -y(t)*abs(y(t)); 
ic1:= y(0) = 1.0, D(y)(0) = 0.; 
dsol1:= dsolve({deq1, ic1}, numeric, range= 0 .. 10); 

        j-> plots:-odeplot(
            [t, eval[recurse](diff(y(t), [t$j]), {deq1})],
            labels= [t, ``],
            title= typeset(diff(y(t), [t$j]))


In your first piecewise command, you've misspelled eps as esp. In your second piecewise, you've done that again, and you've misspelled beta as betta.

In  your last command, t_atan2 means nothing to Maple, unless you've defined it elsewhere. You should replace it with arctan. Maple's arctan can process 1 or 2 arguments.

Using the command subs as you've shown doesn't "set" anything; it can only do so if its final argument is some expression that contains beta.

I recommend against using and (or other logic operators) in the conditions of piecewise. The processing of the conditions proceeds left to right, thus each subsequent condition automatically includes the negation of the disjunction of all previous conditions. Making use of this fact, you almost never need to use and, etc. In other words, 
piecewise(cond1, expr1, cond2, expr2, cond3, expr3)
is logically equivalent to
    cond1,  expr1,
    not cond1 and cond2,  expr2,
    not (cond1 or cond2) and cond3,  expr3,
    not (cond1 or cond2 or cond3),  0

So, I recommend this:
W:= beta-> piecewise(
        beta < Pi/2-theta__M, 
        beta <= Pi/2+theta__M, 
        beta <= Pi-theta__M, 
W(arctan(diff(pho(x,y),y), diff(pho(x,y),x)));


What you're calling a "comma-separated vector" Maple simply calls a list. Your x and f are lists. There is no need for you to use Matrix or Vector for what you show above; it can all be done with lists:

x:= [x__1, x__2, x__3]:
f:= x__1: #no need for f to be a list
Gradf:= diff~(f, x);
Hessf:= map2(diff~, Gradf, x); 

What you call Gradf2 is usually called the Hessian of f.

If you still want to use matrices or vectors, a matrix or vector v can be converted to a list by [seq](v).

Don't use evalf. Its sole purpose is to make decimal approximations. Here are a few ways to do what you want:

  1. exp(2*Pi*I*i/3) $ i= 0..2
  2. seq(exp(2*Pi*I*i/3), i= 0..2)
  3. exp~(2*Pi*I/3*~[$0..2])

As we have discussed here several times in the past, you are not using mod correctly. Use irem instead of mod.

You have found a bug in Maple, and I'll give you a workaround in a moment. But your code itself also has 1 serious bug and 1 serious shortcoming, and it amazes me that these didn't confound the issue of the true bug that you found.

Bug: Change
f:= x-> piecewise(0 <= t, 1, t < 0, t - 1)
f:= t-> piecewise(0 <= t, 1, t < 0, t - 1)

Shortcoming: Change
F:= x-> int(f(t), t= 1..x)
F:= unapply(int(f(t), t= 1..x), [x])

Workaround: Change the unapply to
F:= unapply(int(f(t), t= 1..x, method= FTOC), [x])

So, why use unapply? Because contains a symbolic operation, the integration, that can be performed independently (as well as once-and-for-all) of F's parameter x.


The Taylor or Maclaurin series for a function of a square matrix argument is the same as the series for the same function with a complex argument. The general non-commutativity of matrix multiplication plays no role in this because all powers of a single square matrix do commute. Maple can find the general term of the series for many fairly simple functions and express those series in sigma notation, like this:

convert(exp(A), FormalPowerSeries); #The second argument can be abbreviated FPS.

The commands series or taylor (essentially the same thing) can give you any finite number of terms of the series for nearly any function that can be represented as a power series, but they won't give the general term:

series(exp(A), A=0, 9);

Here is a procedure for computing the nth partial sum of the series for exp(A). The simplicity of this formulation amazes me! It avoids computing high powers of the matrix (which can easily lead to numeric overflow) by combining the updating of the matrix power and the updating of the factorial denominator with the single operation Ak.= A/k. This is an updating assignment operator, introduced in Maple 2018. It's equivalent to Ak:= Ak.(A/k).

#Computes the nth partial sum of the Maclaurin series for exp(A):
Exp:= (A::Matrix(square), n::nonnegint)->
local k, Ak:= rtable(identity, rtable_dims(A), rtable_options(A));
    Ak + add((Ak.= A/k), k= 1..n)
Digits:= 32:
LA:= LinearAlgebra:
Display:= x-> (x, print(evalf[3]~(x)))[1]:
A:= (Display@LA:-RandomMatrix)(
    (6,6), generator= rand(-2.0..2.0), datatype= sfloat
              [ 0.621    1.09     1.55  0.810  -0.464    0.198]
              [                                               ]
              [ -1.41   0.573  -0.0609  -1.66  0.0230    0.816]
              [                                               ]
              [ 0.689  -0.178     1.51  0.384    1.73  -0.0984]
              [                                               ]
              [  1.60  -0.448    0.815   1.79   -1.90    -1.53]
              [                                               ]
              [-0.209  -0.588    -1.19  -1.27   0.338     1.36]
              [                                               ]
              [ -1.73   0.935    0.441  -1.87   0.344  -0.0294]

E1:= Display(CodeTools:-Usage(Exp(A, 99))):
memory used=10.54MiB, alloc change=0 bytes, 
cpu time=62.00ms, real time=68.00ms, gc time=0ns

                [ 3.89     1.54   6.33   2.37   1.29  -0.195]
                [                                           ]
                [-20.6    0.996  -18.2  -23.9   7.80    11.1]
                [                                           ]
                [0.660   -0.262   3.31  -1.55   4.30    2.03]
                [                                           ]
                [ 30.8    0.288   29.1   37.0  -12.3   -17.5]
                [                                           ]
                [-14.5   -0.328  -14.5  -17.3   5.83    8.49]
                [                                           ]
                [-21.0  -0.0767  -18.4  -24.7   8.69    12.1]

E2:= Display(CodeTools:-Usage(LA:-MatrixExponential(A))):
memory used=2.85MiB, alloc change=0 bytes, 
cpu time=16.00ms, real time=14.00ms, gc time=0ns

                [ 3.89     1.54   6.33   2.37   1.29  -0.195]
                [                                           ]
                [-20.6    0.996  -18.2  -23.9   7.80    11.1]
                [                                           ]
                [0.660   -0.262   3.31  -1.55   4.30    2.03]
                [                                           ]
                [ 30.8    0.288   29.1   37.0  -12.3   -17.5]
                [                                           ]
                [-14.5   -0.328  -14.5  -17.3   5.83    8.49]
                [                                           ]
                [-21.0  -0.0767  -18.4  -24.7   8.69    12.1]

#absolute and relative errors:
Display(<[Error__||(abs,rel)]=~ LA:-Norm~([E1-E2, (E1-E2)/~E2])>):
                          [                    -29]
                          [Error__abs = 1.72 10   ]
                          [                       ]
                          [                    -30]
                          [Error__rel = 3.99 10   ]

#My procedure Exp can handle symbolic matrices just as well. Here is MacDude's
#example done by it:
Exp(<a, b; c, d>, 2);
                [        1  2   1           1       1      ]
                [1 + a + - a  + - b c   b + - a b + - b d  ]
                [        2      2           2       2      ]
                [                                          ]
                [     1       1                1       1  2]
                [ c + - c a + - d c    1 + d + - b c + - d ]
                [     2       2                2       2   ]

If the square matrix A is diagonalizable (which is true for most numeric matrices in practice), then there's a better way than partial sums to compute f(A): If A = P^(-1).D.P where is a diagonal matrix, then f(A) = P^(-1).f(D).P where f(D) can be computed simply by applying the complex-valued f to the entries on the diagonal. This is why the MatrixExponential command shown above is faster.

Many unexpected things can happen with alias because it acts at the level of what we call automatic simplification, which is one of Maple's first steps in the evaluation of an expression and often overlooks the intended meaning of an expression that would otherwise be obvious at the deeper levels of evaluation[*1]. Thus, alias must be used with great care. It is often one of the first commands taught to a student just learning Maple, as it was for me. That is perhaps not a good idea. From years of answering Questions here, I've seen that problems often occur when trying to mix alias with indexed structures (either explicitly indexed such as a[1] or indexed-in-essence such as a1a2)[*1].

Here's an alternative to alias that still gives you all the benefits of abbreviation:

n:= 3:
A:= [a||(1..n)](r);
dA:= diff(A,r);

Although it's not an issue in this case, one should also keep in mind that variable names constructed by concatenation are always global, regardless of whether the construction is done by ||catnprintfparse, or convert(..., name), regardless of whether the base name (a in this case) is local or global, and regardless of whether the construction occurs inside a procedure.

[*1] This paragraph from the help page ?alias explains a lot of this:

  • Because aliases are resolved at the time of parsing the original input, they substitute literally, without regard to values of variables and expressions.  For example, alias(a[1]=exp(1)) followed by evalf(a[1]) will replace a[1] with exp(1) to give the result 2.718281828, but evalf(a[i]) will not substitute a[i] for its alias even when i=1.  This is because a[i] does not literally match a[1].

@Zeineb All 4 definite integrals from your most-recent worksheet compute very quickly to exact real values in Maple 2020. In your Maple, try doing it by limits, such as:

I3:= int(f*P3, x); 
simplify(eval(I3, x= 2) - limit(I3, x= 0, right));

This is also worth trying, and is less typing if it works:

int(f*P3, x= 0..2, method= FTOC); 

My guess is that the two methods above are equivalent. FTOC stands for Fundamental Theorem Of Calculus.

Most of the y-values have imaginary parts of significant magnitude.

4 5 6 7 8 9 10 Last Page 6 of 332