Regarding your first question,

> analyticReconstruct:=proc(U,X,Y,Z)
>   simplify(eval(U-I*int(diff(U,Y),X),X=Z-I*Y));
> end proc:

> analyticReconstruct(x^2-y^2,x,y,z);
                                       2
                                      z
 
> analyticReconstruct(exp(x)*cos(y),x,y,z);
                       exp(z - y I) (cos(y) + sin(y) I)
 
> convert(%,exp);
                                    exp(z)

> analyticReconstruct(exp(x)*sin(y),x,y,z);
                      -exp(z - y I) (-sin(y) + cos(y) I)
 
> convert(%,exp);
                                   -I exp(z)
 
> evalc(Re(eval(%,z=x+I*y)));
                                 exp(x) sin(y)

acer

1) You do not have to assign the values in the results from fsolve (to their respective variable names) in order to use them further on. In fact, it's often more convenient to have the variables continue to be unassigned. That's why it is not done by default.

Personally, I think that automatic assignment of solve or fsolve or Optimization results would make the system highly unusable more generally. There's more in this world than just one person's type of problem. A typical PITA would be if one wanted to create other equations or expressions involving those same names, after the first fsolve use. That would necessitate unassigning those names, which is another type of PITA.

> sol:=fsolve({sin(x*y),x+y},{x,y});
                  sol := {x = 26.58680776, y = -26.58680776}
 
> sol;
                      {x = 26.58680776, y = -26.58680776}
 
> eval( 2*x-y/16, sol ); # use the equations in sol
                                  54.83529100
 
> x, y; # handily, x and y are still unassigned
                                     x, y

> neweq := x - tan(y); # easy to get, here
                              neweq := x - tan(y)
 
> assign(sol); # the other way to do it
> x, y;
                           26.58680776, -26.58680776
 
> 2*x-y/16;
                                  54.83529100

> neweq := x - tan(y); # now difficult to get, with x,y, symbols
                             neweq := 35.11465548

> unassign('x','y'):
> neweq := x - tan(y);
                              neweq := x - tan(y)

2) If you lprint the result of the first fsolve call in your worksheet, you can see that the result does actually contain the atomic identifiers. It is not fsolve which converts any variable names to indexed names. It is the GUI that does it, sometimes when you cut and paste.

acer

Generally speaking, by "recursive procedure" is meant a procedure which calls itself.

For example,

> bisection:=proc(f,a,b,eps)
> local c;
> c:=evalf((a+b)/2):
> if abs(evalf(b-a))<=eps then return c; end if;
> if signum(f(c))=signum(f(a)) then
>  bisection(f,c,b,eps);
> else
>  bisection(f,a,c,eps);
> end if;
> end proc:
>
> bisection(sin,1,4,1.0e-7);
                                  3.141592667
 
>
> evalf(Pi);
                                  3.141592654

You can also use procname instead of the name bisection itself, to refer to that current procedure within its own body.

There is lots of room for other improvements in that example. You should keep in mind the danger of non-terminating code.

acer

> ee:=4*cos(theta)^3-3*cos(theta) = cos(3*theta);
                                 3
               ee := 4 cos(theta)  - 3 cos(theta) = cos(3 theta)
 
> simplify(lhs(ee)-rhs(ee),trig);
                                       0
 
> simplify(lhs(ee)-rhs(ee));
                                       0

> expand(ee); evalb(%);
                      3                              3
          4 cos(theta)  - 3 cos(theta) = 4 cos(theta)  - 3 cos(theta)
 
                                     true

acer

Matrix(5, 5, (i,j)->`if`(i=j,f(i,j),0));

Matrix(5,5,scan=diagonal,[seq(f(i,i),i=1..5)]);

LinearAlgebra:-DiagonalMatrix([seq(f(i,i),i=1..5)]);

M:=Matrix(5,5):
ArrayTools:-Copy(5,Vector(5,i->f(i,i)),0,1,M,0,5+1);
M;

acer

The topic might be Markov chains. See here for a note of this in terms of so-called balance equations. Now, most of us are familiar with the question of diagonalizing a Matrix, in the usual linear algebra sense that Patrick showed. But for Markov chains, that same notation is used for qiute a different concept -- related to reversibility of the transition matrix.

There are several sufficient conditions for reversibility of a transition matrix. For example, there are sufficient conditions related to equal entries (excepting the diagonal)  for any given row (or column), but neither of those can hold here, it appears. Also, a transition Matrix A is reversible if there exist symmetric B and diagonal P such that A=P.B.P^(-1), from which it follows I think that any symmetric matrix is reversible. Hence x=1 should provide at least one solution for the first posted example. And a transition matrix cannot be reversible if off-diagonal nonsymmetric entries are zero, so x<>0.

Now some possible good news. For a 3x3 matrix only one (balance) equation needs to be checked, to ascertain reversibility. Check whether A[1,3]*A[3,2]*A[2,1] is equal to A[1,2]*A[2,3]*A[3,1]. Here, that means that 2=2*x for the first posted example, and so x=1 is (also) necesssary, I believe.

I suspect that P^(-1).A.P does not mean "diagonalize A" in the usual linear algebra sense. I suspect that it means find a symmetric B=P^(-1).A.P where P is diagonal. So the second example is fun.

acer

You might create and use a new system of units, as the system SI augmented with the unit 'bytes'. Below, such an augmented system is created and used (for simple examples at least), with the name OPSI.

> restart:

> thing:=3*Unit(gibibytes) + 7*Unit(mebibytes);
                     thing := 3 [gibibyte] + 7 [mebibyte]
 
> combine(thing,'units');
Error, (in Units:-Standard:-+) cannot represent the given unit
in the system `SI`

> simplify(thing);
Error, (in Units:-Standard:-+) cannot represent the given unit
in the system `SI`

> Units:-Standard:-`+`(3*Unit(gibibytes), 7*Unit(mebibytes));
Error, (in Units:-Standard:-+) cannot represent the given unit
in the system `SI`

> Units:-AddSystem(OPSI,Units:-GetSystem('SI'),'bytes');
> Units:-UseSystem(OPSI);

> combine(thing,'units');
                               3228565504 [byte]
 
> simplify(thing);
                               3228565504 [byte]
 
> Units:-Standard:-`+`(3*Unit(gibibytes), 7*Unit(mebibytes));
                               3228565504 [byte]
 
> with(Units:-Standard):
> 3*Unit(gibibytes) + 7*Unit(mebibytes);
                               3228565504 [byte]

acer

A:=Matrix(5,5,(i,j)->1/(1/2+i+j));

acer

One way is to just use the Matrix constructor itself to create the new result.

Mnew := Matrix(5,1,(i,j)->(c^2*s[i,j]^2-x[i,j]^2-y[i,j]^2-z[i,j]^2));

Another way might be to use `zip` or elementwise `^` to get the Matrices with squared entries, and then just do usual addition & subtraction (or scaling by c^2) with those.

Have you considered using Vectors instead of 5x1 Matrices for s, x, y, and z? If you did, then you could use the Vector constructor in a similar way,

Vnew := Vector(5,(i)->(c^2*s[i]^2-x[i]^2-y[i]^2-z[i]^2));

Actually, judging by the error message you showed about, it may even be that you are already using Vectors and not 5x1 matrices as initially claimed. If so, then that Vector constructor call above could serve.

I don't have Maple in front of me, but what would the new elementwise syntax be for this? Does it work ok for caret? c^2*s^~2 - x^~2 - y^~2 -z^~2

acer

If you are calling dsolve,numeric inside a loop or process with many iterations, then memory allocation can be a problem. If that is the difficulty you're seeing, and if your ODEs are IVPs which differ only by some parameters in the formulae then you may be able to pull the dsolve call outside the loop, using its parameters option.

There is a simple example of the difference in performance and resources here. But I've seen a GlobalOptimization example where using the technique brought memory allocation down from many GB to a mere 40MB.

acer

Your expressions appear to have inert Int calls in them. Try value or evalf, not eval.

 (**) expr:=Int(x*M,x=0..1/(3*M));
                                         1
                                        ---
                                        3 M
                                       /
                                      |
                             expr :=  |     x M dx
                                      |
                                     /
                                       0

 (**) eval(expr,M=1);
                                     1/3
                                    /
                                   |
                                   |     x dx
                                   |
                                  /
                                    0

 (**) value(subs(M=1,expr));
                                     1/18

 (**) evalf(subs(M=1,expr));
                                 0.05555555556

acer

I'm sure that someone will provide a good response to your "area" query. I'd like to make a few side comments.

The code you posted for routine MyMandel looks similar to what is in one of the "code edit regions"of the FractalFun worksheet available here on the Application Centre.

That routine is (relatively) very slow when run using Maple's normal interpreter (which is what that worksheet does, using LinearAlgebra:-Map to populate a Matrix to be used to form an image). The whole construction would run considerably faster if edited to utilize Maple's evalhf interpreter to use MyMandel as initializer in a Matrix constructor call. And it would run faster still if it used Maple's Compiler:-Compile facility on the MyMandel routine. Much faster still would be a Compiled version of MyMandel which accepted a Matrix as a new argument (to be acted upon inplace).

It's surprising thing that the version of the demo up on the Application Centre is so slow, and doesn't make use of  Compile. This is precisely the sort of numeric routine for which the Compiler seems designed. (Maple on 32bit MS-Windows is bundled with the Watcom C compiler, which Compile uses. Every other platform can use a free compiler, such as gcc. So availability is not an issue.) And the demo doesn't even benefit from evalhf, which has been built into Maple versions for over 10 years and is a standard way to accelerate one's floating-point Maple routines.

Sure, one might alternatively use Maple's new Threads package to populate the Matrix (image data) for this embarrassingly parallelizable task. But that would get only a few-fold speedup (linear w.r.t. the number of cores, if one were lucky...). Compiling the whole thing, on the other hand, brings 100s of time speedup immediately, with little extra effort. I suspect that this is the sort of reasoning behind Alec's comment here.

Here is a version one might compare with that in the demo worksheet.

restart:

MyMandel := proc(N::integer[4], X::Vector(datatype=float[8]),
                 Y::Vector(datatype=float[8]),
                 M::Array(datatype=float[8]))
local ii::integer[4],jj::integer[4],i::integer[4],new::integer[4],
      zvalr::float[8],zvali::float[8],newzvalr::float[8],
      cvalr::float[8],cvali::float[8];
  for ii from 1 to N do
    for jj from 1 to N do
      zvalr:=X[jj]; zvali:=Y[ii];
      for i from 2 to 25 do
        newzvalr:=zvalr^2-zvali^2+X[jj];
        zvali:=2*zvalr*zvali+Y[ii];
        zvalr:=newzvalr;
        if zvalr^2+zvali^2 > 4 then
          M[ii,jj]:=i; break;
        end if:
      end do:
    end do;
  end do;
  NULL;
end proc:

t1, ba := time(), kernelopts(bytesalloc):

pts := 400:
(X1, X2, Y1, Y2) :=-2.0, 0.7, -1.35, 1.35:

X := Vector(pts, i->X1+(X2-X1)*(i-1)/(pts-1), datatype=float[4]):
X := Vector(X,datatype=float[8]):
Y := Vector(pts, i->(Y1+(Y2-Y1)*(i-1)/(pts-1)), datatype=float[8]):

pMyMandel := Compiler:-Compile(MyMandel):
Val := Array(1..pts, 1..pts, datatype=float[8]):
pMyMandel(pts,X,Y,Val):

time()-t1,kernelopts(bytesalloc)-ba;

img:=ImageTools:-Create(Val,fit=true):
ImageTools:-View(img);

For this particular fractal image's creation, another factor of two speedup could be attained by a cheap symmetry trick. But it's hardly necessary, I think. If anyone wants to post hard timing numbers, to compare with the original...

acer

Is this of some use? (You didn't mention whether variables K, s, i1, i2 and i3 were to be part of the solution or not. It may be that those cannot all be eliminated.)

> eqns:=[Vo = K*(Vb-Vo), i3 = C3*Vb*s, Va-Vb = R3*i3,
>        i2 = C2*s*(Vo-Va), Vi-Va = R1*i1, i3 = i1+i2]:

> Vo=eval(Vo,eliminate(eqns,{Vo,Vb,Va,i1,i2,i3})[1]);

Vo = K Vi/(1 + s R3 C3 + R1 C3 s + K + s R3 C3 K + K R1 C3 s + C2 R1 s
 
           2                2
     + C2 s  R3 C3 R1 + C2 s  R3 C3 K R1)

acer

Aside from other answers, I would suggest that it's generally not good practice to use a protected name such as `sum` as a local variable.

For fun, compare with,

f := proc(a)
local s;
  s := a mod 9;
  if s=0 then 9 else s; end if;
end proc:

acer

It's Matrix(), not MATRIX(), to construct them.

The "end if" and "end do" statements are all mixed up and interleaved. Those have to be sorted into a logically valid layout.

Loops like this simply repeat the exact same command each time through, since the left hand side of the assignment, X_0, always refers to the same thing. It seems very likely that's not what you intend.

for i from 1 to n do
  X_0:=0
end do:

acer