~ was introduced specifically for that, see ?operators,elementwise .

It works in the original example as

ln~(u-v);

                            [    Pi I    ]
                            [            ]
                            [ln(2) + Pi I]
                            [            ]
                            [ln(3) + Pi I]

or, if you want to use f,

f~(u,v);

                            [    Pi I    ]
                            [            ]
                            [ln(2) + Pi I]
                            [            ]
                            [ln(3) + Pi I]

similarly for 3 or more Vectors,

u1,u2,u3:=<1,1,1>,<1,2,3>,<2,3,4>:
ln~(u1+u2)-u3;

                            [ ln(2) - 2 ]
                            [           ]
                            [ ln(3) - 3 ]
                            [           ]
                            [2 ln(2) - 4]

f:=(x,y,z)->ln(x+y)-z:
f~(u1,u2,u3);

                            [ ln(2) - 2 ]
                            [           ]
                            [ ln(3) - 3 ]
                            [           ]
                            [2 ln(2) - 4]

Alec

In Maple 9, you can do

t1b:=map(rhs@op,[t1a]);

In recent Maple versions, the following also works,

t1b:=(rhs@op)~([t1a]);

Alec

For example,

f:=proc(n,p,alpha,c)
local k;
k:=1/(Zeta(alpha)-Zeta(0,alpha,n+1));
fsolve(1-sum(1/(k*p*c+m*i^alpha), i = 1 .. n),m=1,0..1)
end:

f(10,1/2,3,1);

                             0.8637165918

f(1000,1/4,7,1);

                             0.7631056622

You can also use the following command for plotting, to check visually that the solution is unique,

g:=proc(n,p,alpha,c)
local k;
k:=1/(Zeta(alpha)-Zeta(0,alpha,n+1));
plot(1-sum(1/(k*p*c+m*i^alpha), i = 1 .. n),m=0.1..1)
end:

Alec

You could use Wolfram|Alpha, or try Axiom, FriCAS, OpenAxiom, Reduce, Maxima, or Sage online (click on the corresponding link in the Try Online section at the lhs), or Sage in one of its online servers.

Magma calculator is also available online. The input in this example (with 20 digits) is

R<x>:=PolynomialRing(ComplexField(20));
Roots(227.0207672*x^4 - 703.4677734*x^3 + 845.7313843*x^2 - 465.6021423*x + 98.41660309);

Which produces the following output:

[ <0.78168742284054963613 + 0.35390240352655021979*$.1, 1>,
<0.78168742284054963613 - 0.35390240352655021979*$.1, 1>,
<0.79031807704503989295, 1>, <0.74500056719701531879, 1> ]

1s at the end are root multiplicities. The difference with Maple is that Maple produces all 20 digits correct, and Magma calculator gives last 3 digits wrong - by design.

Alec

PS How did you get this equation if you don't have any math software installed? -Alec

PPS About Maple's Online Math Oracles - I couldn't make them working in Firefox and IE in usual mode, but the Solving Equations applet had appeared in IE in compatibility mode. I couldn't paste the equation in the input box there, but out of curiosity, I typed it manually. It was solved fast (under 5 seconds) and the (10-digits) results were correct. However, I wasn't able to copy them from there to post here - obviously, users are expected to write them down if they need them. -Alec

PPPS If nothing else is installed - Python should be installed, and numpy. It's quite simple (and correct) in Python with numpy,

>>> import numpy as np
>>> np.roots([227.0207672, -703.4677734, 845.7313843, -465.6021423, 98.41660309])

array([ 0.78168742+0.3539024j,  0.78168742-0.3539024j,
        0.79031808+0.j       ,  0.74500057+0.j       ])

Alec

Check the view option. For example,

plots:-pointplot([[1,2],[2,3]],view=[0..4,0..4],symbolsize=30);

Alec

Some automation could be done as follows,

numtovec:=L->[seq]([[Re,Im]~(i[..-2])[],color=i[-1],_rest],i=L):
vecplot:=proc({width:=0.1,head_length:=0.5})
plots:-display((plots:-arrow@op)~(numtovec(select(type,[_passed],list),_options)),
scaling=constrained,axes=box,gridlines) end:

Then, the plotting can be done as

vecplot([z1,red],[z1,z2,blue],[z1+z2,yellow]);

or with adding width and/or head length,

vecplot([z1,red],[z1,z2,blue],[z1+z2,yellow],width=0.2,head_length=0.6);

Additional options could be added to the vecplot command, if necessary, next to the width and head_length.

Alec

One can use Joe Riel's NumberToWords module together with the following procedure modifying a number written in words to its ordered form.

ordered:=proc(s)
local T,w;
T:=table(
["one"="first"
, "two"="second"
, "three"="third"
, "five"="fifth"
, "eight"="eighth"
, "nine"="ninth"
, "twelve"="twelfth"
, "twenty"="twentieth" 
, "thirty"="thirtieth"
, "forty"="fortieth"
, "fifty"="fiftieth"
, "sixty"="sixtieth"
, "seventy"="seventieth"
, "eighty"="eightieth"
, "ninety"="ninetieth"]);
w:=StringTools:-Split(s," -")[-1];
if T[w]::string then cat(s[..-length(w)-1],T[w])
else cat(s,"th") fi
end:

This procedure works as

ordered(NumberToWords(23970927409842));

  "twenty-three trillion nine hundred seventy billion nine hundred\
         twenty-seven million four hundred nine thousand eight hun\
        dred forty-second"

Now, the A049525 can be constructed as

A049525:=proc(n::posint)
local s,i,r;
uses StringTools;
r:=Array(1..n);
r[1]:=1;
s:="Iisthefirst";
for i from 2 to n do 
r[i]:=SearchText("i",s,r[i-1]+1..-1)+r[i-1];
s:=cat(s,DeleteSpace(SubstituteAll(ordered(NumberToWords(r[i])),"-","")))
od;
seq(i,i=r)
end:

A049525(100);

  1, 2, 8, 19, 25, 41, 51, 56, 61, 66, 71, 76, 81, 86, 91, 103, 115,

        120, 126, 131, 137, 142, 148, 164, 178, 201, 222, 238, 243,

        259, 307, 323, 351, 367, 405, 410, 432, 446, 451, 494, 510,

        515, 532, 555, 588, 615, 631, 636, 652, 664, 680, 691, 700,

        712, 723, 734, 739, 744, 755, 761, 767, 777, 786, 797, 802,

        807, 818, 823, 828, 838, 849, 859, 870, 880, 884, 889, 899,

        905, 962, 979, 1003, 1008, 1048, 1053, 1070, 1075, 1092, 1142,

        1148, 1165, 1178, 1196, 1215, 1227, 1237, 1256, 1260, 1278,

        1284, 1297

Alec

In Standard Maple, type

coprod

and click Ctrl and Space simultaneously.

Or open Large Operators palette and click on that symbol.

To attach upper and bottom limits, start from the Layout palette. Click on the A with b at the bottom and c at the top - when it appears in the worksheet, A is higlighted, so either click the symbol in the Large Operators palette, or type coprod, then Ctrl+Space, and the symbol will appear. Similarly for limits - delete b and c and type whatever you want there.

Alec

For the latter sum, Maple gives the correct answer,

sum(1/n^(D-2),n=1..infinity) assuming D>3;

                             Zeta(D - 2)

Assuming D::posint is not good enough here because the series diverges for D≤3.

For the integral, the correct answer with assumptions D>2, beta>0 is GAMMA(D-1)*Zeta(D-1)/beta^(D-1) .

For D=2 and beta>0, the integral diverges (at 0) because for small p exp(beta*p)=1+beta*p+... and the integral of 1/p diverges at 0; and similarly for D<2 the integral diverges (at 0).

Alec

PS D in Maple is reserved for function derivatives, so it would be better to use another letter instead of it (but not I, Beta, Pi, Psi, Zeta, or gamma) -Alec

PPS A humorous thing is that Maple is not able to evaluate the integral form of Zeta in its own FunctionAdvisor,

FunctionAdvisor(integral_form,Zeta(z+1))[2];

  [Zeta(z + 1) =

                                   infinity
                                  /                z
                   1             |              _k1
        -----------------------  |          ------------ d_k1,
        /     1  \               |          1 + exp(_k1)
        |1 - ----| GAMMA(z + 1) /
        |      z |                0
        \     2  /

        And(0 < 1 + Re(z))]

value(indets(%,specfunc(anything,Int))[]) assuming z>0;

                       infinity
                      /                z
                     |              _k1
                     |          ------------ d_k1
                     |          1 + exp(_k1)
                    /
                      0

Hovever, this integral can be evaluated using the following,

convert(Zeta(k+1),Int)assuming k>0;

                                    infinity
                                   /                k
                   1              |              _k1
        ------------------------  |          ------------ d_k1
              (-k)                |          1 + exp(_k1)
        (1 - 2    ) GAMMA(k + 1) /
                                   0

isolate(Zeta(k+1)=%,Int);

     infinity
    /                k
   |              _k1
   |          ------------ d_k1 =
   |          1 + exp(_k1)
  /
    0

                          (-k)
        Zeta(k + 1) (1 - 2    ) GAMMA(k + 1)

Now, we can reduce our original integral to this one as follows,

J1:=int(p^k/(exp(p)-1),p=0..infinity):
J2:=int(q^k/(exp(2*q)-1),q=0..infinity):
IntegrationTools:-Change(J1,p=2*q);

                        infinity
                       /          (k + 1)  k
                      |          2        q
                      |          ------------ dq
                      |          exp(2 q) - 1
                     /
                       0

so J2 = J1/2^(k+1). Now

IntegrationTools:-Combine(J1-2*J2);

                         infinity
                        /              k
                       |              p
                       |          ---------- dp
                       |          exp(p) + 1
                      /
                        0

The beta can be added just by changing variables.

Alec

In addition to the way suggested by Markiyan Hirnyk, that could be done also as

degree(f-eval(f,x=0));

                                  6

The usage is approximately the same,

p:=expand(randpoly(x,dense,degree=1000,coeffs=(()->randpoly(y,dense,degree=2000)))+
randpoly(y,dense,degree=1000,coeffs=(()->randpoly(z,dense,degree=2000)))):

CodeTools:-Usage(degree(p-eval(p,x=0)));
memory used=121.58MiB, alloc change=121.60MiB, cpu time=7.50s, real time=7.64s

                                 3000

CodeTools:-Usage(degree(select(has, p, z)));
memory used=218.90MiB, alloc change=128.04MiB, cpu time=7.38s, real time=7.44s

                                 3000

Alec

You could convert your piecewise function to Heaviside form,

convert(piecewise(c(x,t)>v(x,t),c(x,t)*diff(v(x,t),x$2)),Heaviside);

                  / 2         \
                  |d          |
          c(x, t) |--- v(x, t)| Heaviside(c(x, t) - v(x, t))
                  |  2        |
                  \dx         /

Then, with correcting boundary conditions so that c(0,0) and v(0,0) were defined uniquely, that works,

sys:=[diff(c(x,t),t)=v(x,t)*diff(c(x,t),x$2),
    diff(v(x,t),t)=c(x,t)*diff(v(x,t),`$`(x,2))*Heaviside(c(x,t)-v(x,t))]:
bc:={c(0,t)=2,c(1,t)=1,c(x,0)=2,v(0,t)=1,v(1,t)=2,v(x,0)=1}:
sol:=pdsolve(sys,bc,type=numeric,time=t);

  sol := module()
export plot, plot3d, animate, value, settings;
     ...
end module

Alec

For example,

expr:=0.5*x^2+2.35*Pi+sin(0.5*x):
convert(expr,rational);

                         2
                        x     47 Pi
                       ---- + ----- + sin(x/2)
                        2      20

Alec

That can be done with more assumptions,

MultiSeries:-asympt((Q*(R*S+1-S)/(1-R*R3*k*(1-R1)/C/Pec_i*S)*(H_vap+1
/(1-exp(Pec_i)))*exp(C*Pec_i/(R*S+1-S)*(1-R*R3*k*(1-R1)/C/Pec_i*S)*S/k)-H_vap)*(1-exp(Pec_i*
(1-S)/(R*S+1-S)*(1-R*R3*k*(1-R1)/C/Pec_i*S))),Pec_i,1) assuming R>1 , R1>0, C>0 , H_vap
 >0 , H_liq> 0, R3 <=0, 0<S, S<1, k>0, Q>0, Theta0>0, 1-C/(R*S+1-
S)*S/k>0, -C/(R*S+1-S)*S/k+(1-S)/(R*S+1-S)>0;

  /                                           2
  |                           R R3 (-1 + R1) S
  |-(R S + 1 - S) Q H_vap exp(-----------------)
  \                              R S + 1 - S

                                                     \
              (-1 + S) R R3 k (-1 + R1) S        1   |
        exp(- ---------------------------) + O(-----)|
                    (R S + 1 - S) C            Pec_i /

                    /     1 - S            C S      \
                    |- ----------- - ---------------|
                    \  R S + 1 - S   (R S + 1 - S) k/
        /    1     \
        |----------|
        \exp(Pec_i)/

MultiSeries:-asympt((Q*(R*S+1-S)/(1-R*R3*k*(1-R1)/C/Pec_i*S)*(H_vap+1/(1-exp(Pec_i)))*exp(C*Pec_i
/(R*S+1-S)*(1-R*R3*k*(1-R1)/C/Pec_i*S)*S/k)-H_vap)*(1-exp(Pec_i*(1-S)/(R*S+1-S)*(1-R*R3*k*(1-R1)
/C/Pec_i*S))),Pec_i,1) assuming R>1 , R1>0, C>0 , H_vap >0 , H_liq> 0, R3 <=0,
 0<S, S<1, k>0, Q>0, Theta0>0, 1-C/(R*S+1-S)*S/k<0, -C/(R*S+1-S)*S/k+(1-S)/(R*S+1-S)>0;

  /       /                      2                              2
  |       |      R R3 (-1 + R1) S               R R3 (-1 + R1) S
  |-H_vap |Q exp(-----------------) R S + Q exp(-----------------)
  \       \         R S + 1 - S                    R S + 1 - S

                                 2       \
                 R R3 (-1 + R1) S        |
         - Q exp(-----------------) S - 1|
                    R S + 1 - S          /

                                                     \
              (-1 + S) R R3 k (-1 + R1) S        1   |
        exp(- ---------------------------) + O(-----)|
                    (R S + 1 - S) C            Pec_i /

                    /     1 - S   \
                    |- -----------|
                    \  R S + 1 - S/
        /    1     \
        |----------|
        \exp(Pec_i)/

2 other cases give the same answer as the first one.

Now, for the rhs of the 2nd equation,

MultiSeries:-asympt( (exp(Pec_i/(R*S+1-S)*(1-R*R3*k*(1-R1)/C/Pec_i*S)*(1-S))-1)*(H_liq+1
/(1-exp(-C/k*Pec_i/(R*S+1-S)*(1-R*R3*k*(1-R1)/C/Pec_i*S)*S))),Pec_i,1) assuming R>0 ,
 R1>0, C>0 , H_vap >0 , H_liq> 0, R3 <=0, 0<S, S<1, k>0, Q>0,
 Theta0>0, C/(R*S+1-S)*S/k-(1-S)/(R*S+1-S)>0; 
 
        R R3 k (-1 + R1) S (-1 + S)
  exp(- ---------------------------) (H_liq + 1)
              (R S + 1 - S) C

              (-1 + S) Pec_i
        exp(- --------------) + O(H_liq + 1)
               R S + 1 - S

The O-term in the last formula looks strange.

I didn't check the answers manually, so can't tell whether they are correct or wrong.

Alec

If you really meant the squares of the derivatives, then your motion equation's left hand side is the square of diff(h(x1,x2),x1) + diff(h(x1,x2),x2) which works with the first boundary condition as

motion1:=diff(h(x1,x2), x1) + diff(h(x1,x2), x2) = 0:
pdsolve({motion1, ic[1]});

               h(x1, x2) = -ln(2) + ln(2 x1 - 2 x2 - c)

which doesn't satisfy the second condition.

If you meant the second derivatives, it should be written as

motion2:=diff(h(x1,x2), x1,x1) + 2*diff(h(x1,x2), x1,x2) + diff(h(x1,x2), x2,x2) = 0:

which doesn't produce an answer in Maple with given ic.

I didn't get the error you cited in Maple 14.

Alec

The signs are (-1)^k for ceil(exp(sqrt(k))) ≤ n ≤ floor(exp(sqrt(k+1))), with sum(1/n) for these n being asymptotically equal to ln(floor(exp(sqrt(k+1))))-ln(ceil(exp(sqrt(k)))) which is asymptotically equal (ignoring floor and ceil) to sqrt(k+1)-sqrt(k) which is a decreasing function of k with limit 0. 

Alec

Edit. Perhaps, I should mention few more or less obvious and well-known things used in proving of the series convergence above.

First, if we have a series with real terms which are not necessarily positive, something like a₁+a₂-a₃-a₄-a₅+a₆+..., with positive a_k, we can construct an alternating series b₁-b₂+b₃-b₄+... with positive b_k being the sum of a_i in the groups of neighbors in the first series having the same sign, b₁=a₁+a₂, b₂=a₃+a₄+a₅, b₃=a₆+... in this example.

If this alternating series is convergent, then the original series converges, too, because its partial sums S_k with k inside of a group of terms with the same sign are "squeezed" by the partial sums S_m and S_n with m being the last index in the previous group and n being the last index in this group. These S_m and S_n are also partial sums of the alternating series, so they converge if this alternating series converges. Using the Squeeze Theorem, we see that partial sums S_k converge, too (to the same limit), so the original series converges.

Second, replacing series terms during a series convergence proof with their asymptotics, one has to check that the sum of differences between their terms converges.

We did such a replacement twice. First, replacing harmonic sums with logarithms - the convergence of differences between them is well known (and is a part of the usual definition of Euler's constant gamma). Second - getting rid of floor and ceil - in which case the difference between ln(ceil(exp(sqrt(k)))) and ln(exp(sqrt(k))) is not greater than the difference between ln(exp(sqrt(k))+1) and ln(exp(sqrt(k))) which is less than exp(-sqrt(k)) , which follows from the well known series for the ln after combining these 2 logarithms in one. Similarly for ln(floor(exp(sqrt(k+1))) - the (absolute value of the) difference is less than exp(-sqrt(k)). The sum of exp(-sqrt(k)) converges - one can see that from the integral test, for instance, so the sum of the differences converges absolutely by the comparison test.

Third, the series sum((-1)^k*(sqrt(k+1)-sqrt(k)),k=0..infinity) converges by Leibniz rule.

Alec

PS By the way, Maple is not able to find the last sum, equal to (2-4*sqrt(2))*Zeta(-1/2) = (sqrt(2)-1/2)*Zeta(3/2)/Pi, symbolically, but can do that numerically,

evalf(sum((-1)^k*(sqrt(k+1)-sqrt(k)),k=0..infinity));

                             0.7602096252

evalf((2-4*sqrt(2))*Zeta(-1/2),12);

                            0.760209625215

evalf((sqrt(2)-1/2)*Zeta(3/2)/Pi,12);

                            0.760209625217

Alec