See the ?evalf help-page.

acer

You may be able to write a Compilable proc that does what falls inside your i-loop, and possibly also to handle the sorting that follows it.

See here for an example of explicit code to do Statistics:-StandardDeviation, which may be used with Compiler:-Compile.

Such a Compilable proc might also accept a datatype=float[8] Vector parameter SH and act on that inplace, by populating it with the looped results. Of course, the datatype=float[8] Matrix/Array StockReturn would be one of its input parameters.

You could either loop and call such a proc, passing in each permutation, of even better you could pass in the full `per` permutation (as acceptable data for the Compiler) and loop inside the very proc itself.

acer

I've now got a routine that induces barrel, pincushion, and moustache distortion, which seems to work well as a prototype. Right now, it takes about .2-.25 sec on an 800x600 jpg (as a three-layer float[8] Array, in Maple).

I hope to hammer on it a bit more, to attain even less image degradation by using ArrayInterpolation. And then get the inverted operations. And then post the module code.

acer

If V is in fact already a Maple list of (variable) names then you could try plex(op(V)) as argument to Grobner:-Basis.

acer

I'm pretty sure this could be done fine in Maple.

I'd probably start with a simple design like this, beginning with correcting barrel distortion. Write a procedure that accepts as parameters both the input float[8] Array and the output-container float[8] Array (to hold the result), and any optional scalar parameters that characterize the distortion. Then implement the formulae available in several places on the web, to write corrected data into the output-container. Use variable and parameter typing, so that the procedure could be hit with the Compiler. This procedure would return NULL, as it would act inplace on the reuable output-container, minimizing garbage production.

Once that first proc were satisfactory, create a module package with that proc as one export. Write other exports for other distortions.

acer

Use ':-y' for that. It's a reference to the global name y, not the local parameter name (and protected by single uneval quotes, just in case :-y has been assigned).

acer

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